System and method for three-dimensional image rendering and analysis

ABSTRACT

The present invention relates to methods and systems for conducting three-dimensional image analysis and diagnosis and possible treatment relating thereto. The invention includes methods of handling signals containing information (data) relating to three-dimensional representation of objects scanned by a scanning medium. The invention also includes methods of making and analyzing volumetric measurements and changes in volumetric measurements which can be used for the purpose of diagnosis and treatment.

The present application claims priority to two (2) provisionalapplications identified as follows: U.S. Application Ser. No.60/196,208, filed Apr. 11, 2000; and U.S. Application Ser. No.60/253,974, filed Nov. 29, 2000. Each of these earlier filed provisionalapplications are incorporated herein by reference.

The invention was made with government support under R01CA63393 andRO1CA78905 by the National Cancer Institute. The government has certainrights to the invention.

BACKGROUND OF THE INVENTION

The present invention relates to the art of diagnostic imaging, and, inparticular, to methods and systems for conducting highly enhancedthree-dimensional reproduction and analysis of three-dimensional imagesand diagnosis and treatment related thereto.

It is known in the art of image diagnostic systems and methods to scanobjects found in a host patient, retrieve signals resulting from suchscanning, and process the signals to obtain information for possiblediagnosis and treatment. To date, much of the effort has been directedto treatment of such data in two dimensions, and analysis has to acertain extent been constrained by two-dimensional limitations.Consequently, professionals relying on such information have, similarly,been somewhat constrained in their ability to create diagnostic andpossible treatment models.

As a result of the inventors herein having the perspicacity to think andanalyze based on three-dimensional vision, new methods and correspondingsystems are now available for use in various applications. In one of themore basic applications, detection and analysis of objects in the body,especially growths, can be produced and analyzed to a degree of accuracypreviously thought unknown. This is particularly useful especially inthe area of life threatening pathologies, e.g., cancer. One particularlypernicious pathology is lung cancer. The ability to be able to detectearly stage lung cancer has provided strong motivation to improvediagnostic capabilities.

For example, the small pulmonary nodule is a common radiographic findingthat poses a diagnostic and clinical challenge in contemporary medicine.Approximately 170,000 new lung cancers are detected each year in theUnited States alone. Pulmonary nodules are the major indicator of thesemalignant lung cancers, but may also be signs of a variety of benignconditions. Of the all the new lung cancers currently detected each yearin the United States, approximately 90% are fatal, responsible fornearly 160,000 deaths in 1999 alone [46]. In fact, lung cancer is thenumber one cause of cancer mortality, responsible for more deaths in menthan prostate and colon cancer combined, and in women, for approximatelyas many as breast and colon cancer combined [46].

One contributing factor for this low survival rate is the historicallack of early detection, or the detection of the disease at a stagewhere it is most likely to be curable. Until now, the majority of lungcancers are either “symptom-detected,” or found on chest x-rays (CXR) ofasymptomatic patients. In most cases, the tumor is detected at anadvanced stage where opportunity for successful intervention is limited.It has now been shown, through the work of the Cornell Early Lung CancerAction Project (ELCAP) that it may be possible to detect a significantnumber of malignancies at an early stage, where cure rates are much morepromising, through the use of lung cancer screening [28]. Thecornerstone of the ELCAP study is the use of computed tomography (CT)for the detection of small pulmonary nodules and the use of early-repeatCT (ERCT) for evaluation of these nodules, by high-resolution computedtomography (HRCT). It has been reported that low-dose CT has thepotential to detect tumors 4-6 times as frequently as CXR [28]. It istheorized that this may engender a “stage-shift” in tumors detected byCT versus CXR. This additional time for potential intervention holds thepromise of increasing survival from the current 12% to over 70% S-yearsurvival following resection of a stage I cancer [28, 76].

The tools developed in this work concern the application of computervision and image analytic techniques to provide quantitative informationdescribing the characteristics of pulmonary nodules seen in helical CTstudies. They are an effort to capture in a reproducible, mathematicalform the somewhat subjective characteristics traditionally used byradiologists and other physicians to describe these lesions in order todetermine their likelihood of malignancy. Advances in helical CTtechnology are now allowing for the detection of a significant number ofpulmonary nodules less than 1 cm in diameter. At this point, traditionalmethods for nodule characterization (many based on hand-measurement andsubjective evaluation) become less effective, given the small size ofthese lesions. Through the use of computer-aided diagnostic techniques,we have a unique opportunity to assist the radiologist in theincreasingly difficult task of diagnosing these small nodules.Furthermore, as the expected widespread acceptance of CT screening forearly detection of lung cancer increases, computer-aided tools shouldplay a prominent role in the early detection, characterization, and cureof lung cancer in many thousands of patients around the world.

Pulmonary nodules, while most important when they represent malignancy,are frequently the manifestations of a variety of benign conditions suchas hamartoma, benign adenoma, or rarely, sarcoidosis, amyloidosis orinfection. Fungal pulmonary infections (histoplasmosis, cryptococcosis,coccidiomycosis), mycobacterial (tuberculosis) and other bacterialinfections may also, infrequently mimic neoplasms [69].

Establishing a specific diagnosis is important, since surgery is usuallynot needed in the treatment of benign small pulmonary nodules (SPNs). Onthe other hand, SPNs may represent malignancy (cancer) of several types(squamous-cell, adenocarcinoma, large-cell, small-cell,bronchioloalveolar carcinoma, etc.) where expeditious removal bythoracotomy is usually indicated. For example, approximately 5% ofpatients with “small-cell” lung cancer are diagnosed after a chestradiograph (x-ray) shows an SPN. Currently, the definitive diagnosis ismade after the tumor has been removed surgically. If the diagnosis ismade before the operation, clinical decision making is facilitated sincemultiorgan scanning and inediastinoscopy to detect spread or metastasisare performed prior to thoracotomy [85]. In addition, very smallpulmonary nodules may not be detectable in the ordinary chest x-ray orconventional CT and are visible in helical (spiral) CT examinations.Therefore, improvement in the available techniques for the detection(e.g. differentiation from vessels) and differentiating malignant frombenign SPNs are needed and will result in clinical and economicbenefits.

Small pulmonary nodules are usually identified on routine chestradiographs (CXR) in asymptomatic patients or on chest radiographsobtained for other indications. Most nodules have a benign cause, butbronchogenic carcinoma, metastases, and other malignant processes arecommon causes and must be excluded in the differential diagnosis. Asthey are potentially malignant, SPNs require expeditious evaluation.Techniques including computerized tomography (CT) have been used for thedifferential diagnosis of pulmonary nodules. Radiologic evaluationincludes chest x-ray (CXR) or computed tomography (CT) for the detectionof nodules, high-resolution CT (HRCT) for further characterization, andCT and magnetic resonance imaging (MRI) to evaluate the mediastinum, andextrathoracic structures.

The radiologic characteristics that differentiate benign from malignantpulmonary lesions include size, the presence of calcification in a“benign pattern,” the shape and surface characteristics, rate of changein size, clinical characteristics of the patient such as age and smokinghistory, etc. Meticulous assessment of the margin of the nodule isnecessary. Useful characteristics for the differentiation of benign frommalignant include, sharpness of margin, the presence of spicules, lengthof spicules, spicules extending to the visceral pleura, pleural tailsign and circumscribed pleural thickening. Even with current criteriahowever, the best sensitivity and specificity that can be achieved are85% and 78%, respectively. Further improvements in diagnostic accuracyare needed since at the present time many patients with ambiguousfindings are referred for early surgery [72].

Mediastinal CT is the preferred modality for examining the mediastinumin these patients while magnetic resonance imaging is used selectively,e.g. in patients with superior sulcus tumors who are candidates forsurgery. The differential diagnosis is an important consideration,particularly in patients with small pulmonary nodules, due to theincreased sensitivity in lesion detection, increased specificity andlesion (tissue) characterization (usually done with MILI) [29].Metabolic imaging by positron-emission tomography (PET) provides manyadvantages over conventional anatomically-based cross-sectional imaging,although its use is generally resolution-limited to nodules larger than1 cm in diameter. For routine use in oncology, a detailed assessment ofspecific efficiency of PET is indicated [67]. Mediastinoscopy is doneusing CT for guidance [61]. Fine-needle aspiration biopsy is a usefuldiagnostic tool [95]. Thoracic biopsy is a safe and accurate minimallyinvasive surgical approach to resectable peripheral SPN. Many patientsrequire surgical biopsy especially those with adequate physiologicreserve (that is those who would be able to tolerate surgery) [53].Also, video-assisted thoracic surgery (VATS) often allows one toaccomplish the same goal as the comparable open procedure but with lessmorbidity and a shorter hospital stay. With continued development ofinstrumentation, increasingly complex procedures continue to beaccomplished with this technique [37].

Currently the diagnosis and management of SPNs is based on the basicprinciple that every nodule must be regarded as potentially malignantuntil proven otherwise. Malignant nodules are surgically removed unlessthe patient is in such bad shape from other diseases that he or she cannot survive the surgery or there is proof that the cancer (malignanttumor) originated elsewhere in the body and that the SPN is ametastasis. On the other hand, resection of a benign nodule carries asmall surgical cost. The major benefit derived from resection of thebenign SPNs is that the diagnosis is made by pathological examinationand malignancy is excluded. In addition, once a small pulmonary nodulehas been detected, the diagnosis should be made quickly to avoidspreading of the malignancy.

When a malignant cause cannot be ruled out, the patient's age, smokinghistory, and nodule size must be considered. For moderate and high-riskpatients, an immediate and more invasive work-up is indicated [17].Observation by serial radiographs or CT studies may be the appropriatecourse for patients who are at low risk for malignancy. When theprobability of malignancy is low, the patient may be advised to haverepeated radiologic examinations in order to determine whether the SPNis changing over time before the decision about whether to operate ismade. Therefore, making the diagnosis without surgery either at initialexamination or after consideration of repeated examinations would be ofbenefit [47].

In addition to primary pulmonary pathology, pulmonary nodules may be dueto metastatic disease from cancers of the breast, stomach, pancreas,kidney, bladder, and the genitourinary tract. Computed tomography,especially helical CT, is probably the most sensitive imaging techniquein the identification of pulmonary metastases, because it detects ahigher number of pulmonary nodules compared with other techniques. In apopulation of patients presenting with multiple, solitary pulmonarynodules, or the absence of nodules, helical CT shows 30 to 40% ofsupplementary nodules when compared to conventional scanning techniques.Helical, or spiral, scanning is a relatively new computed tomographictechnique based on the continuous acquisition of volumetric CT dataduring continuous x-ray beam rotation and continuous patient translationthrough the scanner at constant velocity. It has many advantages overconventional CT including rapid image acquisition during asingle-breath-hold scan of the lung, and the ability to obtain axialimage reconstruct ions at arbitrary and overlapping intervals, thusallowing the detection of small lesions that otherwise would beinconspicuous because of respiratory motion or slice averaging. Thisleads to better identification of small pulmonary nodules and to highquality multiplanar reconstructions that can be useful in the study ofmediastinal lymph nodes and the vascular and tracheobronchial spread oflung cancer [10].

In addition to its use in standard radiographic studies, helical CT alsoenables the optimal use of contrast products in which the iodine ispoorly concentrated to study the pulmonary vessels. The continuity, bothanatomically and for the lesions obtained by spiral CT is such that itis now possible to apply to thoracic pathology techniques of multiplanarand three-dimensional reconstruction. If all the information iscontained in conventional transverse imaging is shown slice by slice,not everything is perceived by the observer because the information isinconveniently presented, by deconstructing anatomical and lesionalpicture. In helical CT, reconstructions of the volume may be inspectedfrom a coronal, sagittal, oblique or three-dimensional perspective,furnishing supplementary information to axial images alone [66].

Traditionally, measurements of pulmonary nodules have been made by theradiologist using calipers or a ruler on a chest radiograph. The goalwas to make an approximate measure of the nodule size by estimating thediameter seen in the projection captured on the chest film. The sizescould be recorded, along with subsequent measures to establish whetherthe nodule was growing and at what rate. This practice has carried overinto the reading of many CT examinations. Although computer-basedworkstations are sometimes available for the review of CT studies, thesequential images are often output in “hard copy,” film format andviewed on a light box. Here too, caliper measurements are the norm, withsize being estimated from the image with the greatest cross-sectionalarea.

As a step toward precision and reproducibility of measurement, manydigital review workstations provide “digital calipers,” that allow theradiologist to make measurements using a mouse. While a generalimprovement of caliper-based measurements on film, it should be notedthat the estimate of nodule size is still observer-dependent, for thefollowing reasons: (a) selection of an appropriate cross-section(maximal area) for measurement is variable; (b) nodule cross-sectionsmay not be circular, prompting multiple measurements and averagingschemes; (c) the margin of many nodules is somewhat indeterminate; and(d) window and level settings can markedly change the apparent size of anodule.

These factors suggest that a more automated assessment of nodule sizewould be beneficial to the unbiased, reproducible estimation of size andgrowth. A potentially more important point is that the measurements madeby the radiologist are nearly always two-dimensional in nature. Withcomputer-aided methods, the entire nodule volume may be measured withgreat precision, eliminating the problems of slice selection,non-circularity, and window and level distortion.

Computer vision is the study of techniques for the extraction ofinformation from images. The images may be two-, three-, orhigher-dimensional in nature, and are acquired by varying types of imagesensors, including conventional 2D cameras, video microscopy systems,x-ray film digitization, and computed tomography systems. The underlyingphysics of acquisition also varies, including the absorption of visiblelight, infrared, x-ray, and the sensing of magnetic resonance. Each ofthese may be accomplished in a variety of geometries and projections.

Computer vision is closely allied with its sister disciplines: imageprocessing and computer graphics. Image processing is the application oftechniques for the improvement or modification of an image, generally sothat it may be better perceived by a human observer. Computer graphicsis concerned with basically the inverse problem of computer vision:generation of images based on an informational description.

Most of the techniques now popular in computer vision have beendeveloped for industrial and military applications. Examples includeproduct quality control, and detection of tactical targets in radar,sonar, and conventional images. More recently, computer vision has beenapplied to problems in biology and medicine. Automated analysis of cellmorphology from microscopic images, analysis of anatomical structuresradiographic images, and study of function using advanced imagingmodalities are fast becoming key applications for computer visionmethods.

One notable difference between computer vision in medicine and biologyand traditional computer vision is that in industrial applications, thestructure of objects being studied is frequently known a priori, andalso often has a carefully described geometry. Biological structures,however, while they normally have a particular structure, exhibit a widerange of variation, especially in pathological cases. Even normalanatomical structures defy simple geometric description, as they arethree-dimensional, complex structures that exhibit variation betweensubjects and also often have a dynamic appearance in the same subject.Growth over time, as well as motion during image acquisition make itnon-trivial to obtain easily-modeled data of a particular structure,even in a single subject.

These subtleties motivate techniques that can characterize the size,shape, and other qualities (e.g. density) of a biological structure in aquantitative, reproducible way. This allows for the comparison of databetween subjects, as well as the study of a structure in a singlesubject over time. Such measures are precisely what is needed toeffectively study the progression of a pulmonary nodule from firstdetection to eventual diagnosis, and in many cases, during follow-up aswell.

The purpose of computed tomography (CT) is to display the anatomy of aslice of the body by measuring absorption of x-rays as they pass throughthe body in many trajectories. Imaging usually is done in slicesperpendicular to the axial (head to toe) direction. The resultingmeasurements are manipulated mathematically to obtain two- orthree-dimensional displays.

X-rays are electromagnetic radiation (10⁻¹⁰-10⁻⁸ nm) that are producedwhen an electron beam strikes a metal anode, typically made of tungsten.The anode is heated during the generation of the x-rays and dissipationof heat becomes a limitation in long examinations. Only x-rays in onegeneral direction are allowed to exit the x-ray generator by the use ofa beam collimator made of lead shielding material. As they travelthrough matter, the x-rays are attenuated, i.e., their intensity(related to the number of photons per unit cross sectional area) isdecreased by absorption or scatter. Absorption occurs through thephotoelectric effect, where the x-ray energy is given to an electronthat is freed from the atom. Scatter occurs where x-rays of higherenergy (generated by electrons accelerated by higher voltages, havinghigher frequency) transfer some of their energy to an electron (theenergy necessary to free it from the atom), while the remaining energybecomes a secondary photon (x-ray) that is transmitted in a newdirection. Thus, absorption is useful in diagnostic radiology since itoccurs in a known direction while scatter creates noise due to secondaryx-rays of unknown direction being detected by the x-ray sensor(s).

CT and all other x-ray diagnostic examinations are based on the factthat different structures of the body absorb or scatter x-rays to adifferent degree. The fundamental principle is expressed by Beer's Law

$\begin{matrix}{\frac{\mathbb{d}I}{\mathbb{d}L} = {{- \mu}\; I}} & (1.1)\end{matrix}$describing the attenuation of intensity of electromagnetic radiation asit passes through a homogeneous medium. In this equation, I is theintensity of the radiation (number of x-rays, photons, per unit surfacearea); L is the length of the pathway, and μ, is the linear attenuationcoefficient. This equation merely describes that the x-ray attenuationas it goes through a unit length of the medium is proportional to theintensity of the incident radiation or that each x-ray (photon) hasequal probability of being absorbed while traveling through a giventhickness of a given material as any other x-ray (photon). Thisdifferential equation leads to the following exponential relationdescribing the intensity I of radiation as it passes through a medium:I=I₀e^(−μL)  (1.2)where I is the transmitted intensity through a thickness L. This isderived easily from the fundamental definition of e

$\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}x}( {\mathbb{e}}^{x} )} = {\mathbb{e}}^{x}} & (1.3)\end{matrix}$and the associated chain-rule for differentiable functions of x

$\begin{matrix}{{{\frac{\mathbb{d}}{\mathbb{d}x}( {\mathbb{e}}^{u} )} = {{\mathbb{e}}^{u}\frac{\mathbb{d}u}{\mathbb{d}x}}}{{Therefore},}} & (1.4) \\{\frac{\mathbb{d}I}{\mathbb{d}L} = {\frac{\mathbb{d}( {I_{0}{\mathbb{e}}^{{- \mu}\; L}} )}{\mathbb{d}L} = {{- \mu}\;{I_{0}( {\mathbb{e}}^{{- \mu}\; L} )}}}} & (1.5)\end{matrix}$

FIG. 1 illustrates Equation 1.2, where the intensity of the radiationpassing through the material, I, is a function of the intensity of theincident radiation, I₀, the length of the material, L, and the linearattenuation coefficient of the material, μ. This coefficient isdetermined by the density of the material, the atomic number, and thewavelength (or energy spectrum for polychromatic x-rays) of the x-rayradiation under consideration. Thus, knowing the linear attenuationcoefficient in each voxel of the body would provide information on thedensity and atomic number of the elements at that location.

The process of creating images from x-ray absorption data along avariety of angles and translations is called image reconstruction,involving the solution of a set of equations relating the geometry ofabsorption data seen in each projection. A given CT projection is theresult of x-ray traversal through objects of varying attenuationcoefficients at varying distances. This relationship can be generalizedasI _(n) =I ₀ e−μ ₁ L ₁ ·I ₀ e ^(−μ) ² ^(L) ² . . . I ₀ e−μ _(n) L_(n)  (1.6)where I (I_(n)) is the intensity of the radiation emerging from the bodyand detected by the detector, μ₁, μ₂, . . . μ₃, are attenuationcoefficients of the individual structures traversed by the x-ray and L₁,L₂, . . . L₃, are corresponding distances traveled by the x-rays(lengths). Assuming equal lengths from x-ray source to detectors inhelical CT (L, chosen at the time of image reconstruction), Equation 1.6can be rewritten asI_(n)=I₀ e ⁻ ^((μ) ¹ ^(+μ) ² ^(+ . . . +μ) ^(n) ⁾ L  (1.7)FIG. 2 illustrates the physical system in which Equation 1.7 applies. Inthis case, the incident x-ray beam passes through several materials ofvarying length and attenuation characteristics. The computation of theindividual attenuation coefficients (μ₁, μ₂, . . . , μ_(n)) requires thesolution of an equal number of equations obtained from the differentprojections acquired in the CT procedure.

The usual method for solution of the set of equations and the generationof the CT image is called filtered back-projection, based on the Radontransform. In this method, a uniform value of attenuation is projectedover the path of each ray such that the total attenuation isproportional to the measured attenuation [14, 15, 6]. These values arestored as elements of one vector of the reconstructed matrix. Theprocedure is repeated for each ray sum of the CT while corrections aremade for voxels traversed by the x-ray beam in an oblique fashion. Theassumption that x-ray attenuation is uniform throughout the path of eachx-ray beam (obviously an oversimplification) results in blurring of theresulting image. In spite of this limitation, filtered back-propagationis popular because of its inherent parallelism. It allows processing ofthe data obtained by each ray while data acquisition continues in otherprojections, dramatically improving computational efficiency. Theblurring is decreased by increasing the number of projections and byconvolution with a filter function (convolution kernel). The filteringfunctions depend on x-ray tube geometry, detectors, intended effects(sharpening edges thus enhancing spatial resolution at the expense ofdecreased density resolution versus more refined density resolutionwhile sacrificing spatial sharpness, etc.).

The relative pixel values obtained by this image reconstruction processare not the absolute values of the linear attenuation coefficients ofindividual voxels of the body (μ₁+μ₂+ . . . μ_(n)). The CT numbers,called Hounsfeld Units (HU) in honor of Godfrey Hounsfield (the inventorof CT scanning [30, 31]), are related to the linear attenuationcoefficients as follows:

$\begin{matrix}{{{CT}\mspace{20mu}{number}} = \frac{K( {\mu - \mu_{\omega}} )}{\mu_{\omega}}} & (1.8)\end{matrix}$where μ_(ω), is the attenuation coefficient of water, μ the attenuationcoefficient of the pixel in question, and K a constant. The value of Kmust be large enough (e.g. 200-2000) to accommodate the accuracy of thescanner. For example, consider a CT scanner with a density resolution of±0.5%, or ±1 in 200. In this case, a value of K=200 would besufficiently high to encode density values, as they are typicallyrecorded with integer precision. A larger value of K would expand thescale beyond the accuracy of the scanner.

Helical CT is a newer method of data acquisition in computerizedtomography [83, 40]. Data acquisition is continuous during rotation ofthe x-ray source and the detector (mounted in a toroidal gantry) aroundthe patient who is simultaneously moved in the axial (head to toe)direction through the rotating gantry at precise speeds from 1 mm tomore than 10 mm per second. The resulting helical projections are usedto reconstruct an anisotropic Cartesian space representing the x-rayabsorption at each voxel in the scanned patient volume.

In commonly used instruments, the x-ray beam is collimated to a fan thatdefines the image plane at any given time, and the array of detectorstravels in a circular path around (typically 360°) the patient. Thisapproach provides the opportunity to obtain many projections (andtherefore images) in a short period of time (e.g. in one breath-hold),and the collection of data in a continuous volumetric manner, ratherthan in slices. This allows reconstruction in any position as well asimproved two-dimensional and three-dimensional views. The majorpractical advantage of spiral CT is the ability to cover completeanatomic regions in a single exposure [38]. Isotropic imaging, wherespatial resolution is the same in all three directions, will eventuallybecome a reality [39]. Helical CT has been found useful in the study ofpulmonary nodules where optimally centered studies have been shown toimprove accuracy [40, 88].

The most important parameters of a CT protocol relate to the imagegeometry, reconstruction resolution, and radiation dose. Pitch isdefined as the ratio of the table speed in mm/revolution (assuming a360° reconstruction algorithm) to the collimation width of the x-raybeam. For a given beam collimation, the table may be advanced at a speedthat would move the table through one collimation width per revolution,or 1:1 pitch, or at a higher pitch, reducing the effective radiationdose and scan time at the cost of image quality. For example, a protocolwith 10 mm collimation and 20 nm/revolution table speed would acquireimages with a 2:1 pitch.

The reconstruction resolution is described by several parameters. Thein-plane resolution (the image plane is perpendicular to the scanneraxis), which is equal in the x and y dimensions, gives the resolution ofa pixel in a single 2D reconstructed image. The axial resolution, orresolution in the z dimension, is determined by the slice spacing, orthe physical distance between reconstructed 2D images. This spacing is afunction of the scan collimation, pitch, and reconstruction algorithm.Furthermore, images may be reconstructed at overlapping intervals (e.g.10 mm overlapping slices at 5 mm intervals). The axial resolution ofconventional CT protocols is often 5-20 times more coarse than thein-plane resolution. This results in 3D voxels that are anisotropic (notall dimensions exhibit equal spacing).

Radiation dose is a function of x-ray tube current, typically expressedin mA, tube voltage, expressed in kV (peak measurements are given unitsof kVp), and scan duration. The tube current and voltage directly affectthe quality of reconstructed images, as the signal-to-noise ratio seenat the x-ray detectors is proportional to the dose used. Improvements indetector sensitivity and filtered image reconstruction algorithms arecontinually reducing the minimum dose required to achieve acceptableimage quality at a given scan geometry.

There are three types of CT examinations commonly performed in the studyof pulmonary nodules. The first is a low-dose, low (axial) resolutionscreening scan, aimed at the detection of possible pulmonary nodules inthe lungs. The second, following detection of a pulmonary nodule, is astandard-dose scan of the entire lung volume. The third (and primaryfocus of this work) is a standard-dose, high-resolution scan of theregion containing the detected pulmonary nodule, used for analysis andcharacterization of the nodule itself. A description of the parameterstypically used in the different scan protocols is shown in FIG. 3.

Current CT scanners impose a tradeoff between image quality, slicethickness, and x-ray dose. Many scanners have a gantry rotation speed ofapproximately 1 revolution/second. This implies that a breath hold onthe order of 30 seconds to 1 minute is required to obtain a whole lungscan. More recent scanners, using multislice technology, have reducedthis time to just a few seconds, helping to eliminate respiratory andother motion artifacts. It may be anticipated that future developmentsin CT detector and processing technology may make possible much morerapid and higher resolution scans such that future screening scans maybe taken at high (≈0.1 mm) axial resolution. However, current CADsystems are designed for the constraints of today's scanners whichsuggest a two-resolution approach: a low-dose initial whole-lung scanfollowed (in some protocols) by a high-resolution focused scan ofdetected nodules.

Much research has been done relating to the manual detection andcharacterization of pulmonary nodules by radiologists. In the past 10-15years, work in computer-assisted and automated detection has appeared inthe literature. Until quite recently, the dominant imaging modality forthe study of cancer in the lung has been the chest radiograph, or chestx-ray (CXR). With the increasing use of CT for the detection of theselesions, however, we are now seeing an increase in the use of automatedmethods for the detection of pulmonary nodules in thoracic CT studies.Some work has also been reported dealing with characterization ofpulmonary nodules, but it is only quite recently that advances have beenmade in this area. As characterization of pulmonary nodules from CTscans is the focus of this work, the following sections provide a viewof pertinent related issues in the study of nodules using both CXR, andCT, in the screening setting, and for nodule characterization.

The goal of lung cancer screening is to detect pulmonary nodules,establish their probability of malignancy at an early stage, and therebyimprove likelihood of cure. The traditional detection of pulmonarynodules occurs by their appearance on a routine chest radiograph, or onetaken for unrelated reasons. In the past 10-15 years, research has beendone in developing computer-aided methods for detection of thesenodules.

In the 1970s, the National Cancer Institute (NCI) sponsored three largescreening trials to evaluate the benefits of lung cancer screening usingchest radiography and sputum cytology. These were held at the MayoClinic, Johns Hopkins, and Memorial Sloan-Kettering. [52, 20, 19].Although these studies found that lung cancer screening allowed forearlier detection of lung cancer, better resectability of tumors, andbetter survivorship of the surgery, the overall mortality in thescreened and control groups were similar. Thus, lung cancer screeningwas not recommended as national policy.

Detection of pulmonary nodules has traditionally been done throughvisual inspection of chest radiographs. Since the 1970's, however, muchresearch has been devoted to the automation of this process.Computer-aided detection and analysis systems for pulmonary nodules fallinto several categories. Image processing techniques have been used toimprove contrast and otherwise enhance candidate regions for the visualdetection of nodules [43, 44, 77]. Systems have been developed for theautomated detection of candidate nodules [50, 96, 93, 8, 59]. Some workhas also been done on automated analysis of nodules in CXR images forthe prediction of malignancy [71, 58, 87].

The preprocessing stage of several CAD systems for the detection ofnodules in CXR involves histogram manipulation. Techniques have beenexplored using histogram equalization and high-frequency enhancement[59]. Several studies have described the use of image differencing,which is the subtraction of a “nodule suppressed” image from one that is“nodule enhanced,” to help reduce the influence of other anatomicalstructures appearing in the radiograph [23, 96, 44, 93]. Thisenhancement of CXR images is achieved using one or more matched filtersdesigned to identify the round regions characteristic of many nodules.For example, a disk-shaped filter may be convolved with the image toenhance nodule-like regions. A ring-shaped median filter may be used fornodule suppression. Related work in preliminary segmentation of lungimages has also been done in the areas of lung boundary delineation [91,3] and rib identification [70]. These steps help eliminate regions ofthe image from consideration as potential nodules.

A second stage of histogram manipulation may be performed following theinitial image preprocessing. Multiple thresholding of the preliminaryimages may be performed to identify the incremental contribution of eachintensity band to candidate regions. In this way, suspect regions areiteratively enlarged (by decreasing a threshold) and analyzed [50, 8].Computed tomography (CT) is fast becoming the modality of choice for thedetection and analysis of pulmonary nodules. Studies in lung cancerscreening with CT have been performed in Japan [41, 78] and New York[28].

There have been several prototype systems for performing complete,automated lung CAD [94, 22, 84, 7] (i.e., taking as input whole-lungscans and identifying nodules). Note that these systems involvelow-resolution screening scans and therefore, may be able to detectsmall nodules (<1 cm) but there is insufficient information tocharacterize them. These early attempts employ a multiphase approachinvolving the following three stages: (a) identification of the lungregions in the CT images; (b) separation of candidate nodules fromvessels within the lung regions; and (c) classification of candidatenodules.

For diagnostic classification of small nodules, a focused HRCT scan ofthe region-of-interest is obtained. Prototype CAD systems for thesescans have been developed in which the region-of-interest is specifiedby the radiologist and just the separation and classification stages areautomatic [63, 42].

The lungs may be separated from other structures in the chest by firstperforming some noise reduction through spatial filtering and thenapplying a simple threshold with background removal to the whole-lungscan. As an example, a single slice from a whole-lung screening scan isshown in FIG. 4 and the extracted lung region is shown in FIG. 5. Theentire lung volume can be measured or visualized as shown in FIG. 6. Asimilar method was used by Giger et al. [22], while Brown et al. [71]used an anatomical model with fuzzy set membership based on several softconstraints including threshold, continuity and relative position.Toshioka et al. [84] use geometric filtering to avoid irregularities inthe detected boundary. They also use a second method, basing their lungboundary on predictions from rib locations in order to identify thepresence of very large nodules.

Further analysis of other thoracic structures is also possible. Forexample, the lumen of the trachea and the main bronchi may be traced bya simple seeded region growing algorithm; typical results for screeningdata is shown in FIG. 7. This figure illustrates the limitations of theapproach on screening scans due to their low axial resolution; thebronchial tree and the vascular structure can be extracted to a muchfiner degree with higher-resolution scans. Three-dimensional regiongrowing combined with tree-traversal algorithms has also been used inremoval of the lung vasculature as a preliminary step to noduledetection [90, 16, 79].

Until recently, nearly all of the literature on the characterization ofpulmonary nodules from radiologic images was concerned with theradiographic appearance of the lesions to a trained radiologist. Most ofthe measures of nodule size were made manually using calipers or a ruleron chest radiographs or hardcopy CT images. Measurements of shape anddensity distribution were also described in somewhat qualitative terms.

Some notable work on the characterization and subsequent classificationof pulmonary nodules may be found in [26, 47, 82, 73]. The chiefcharacteristics reported were size, smoothness of edge, presence and(for larger nodules) pattern of classification, spiculation, lobulation,cavitation, and the presence of a pleural tail.

Perhaps more important than shape and density characterization, theestimation of nodule aggressiveness as a function of growth rate hasbeen studied for some forty years. Nodule doubling time, as computedfrom manual measures of nodule diameter, has long been utilizedeffective predictor of nodule malignancy [12, 56].

A variety of techniques have been used to characterize potentialcandidate regions for nodule detection in screening data. Measurement ofobject circularity is often helpful, as this metric may help exclude ribcrossings and other structures [50, 44, 92]. Methods involving gradientanalysis have also been described [50, 49]. These are based on thenotion that the edge gradients of a small nodule-like object shouldprimarily point toward the center of that object, unlike the case withribs and many vessels.

Morphological operations have been explored for the detection of nodulesin CXR [96]. These nonlinear filters impose geometric constraints on thecandidate regions, allowing, for example, round objects (e.g. nodules)to be distinguished from longer, flatter ones (e.g. ribs, heart, etc.).Correlation-based template matching has also been used to identifycandidate nodules, an idea related to matched filter enhancement [8,59].

The selection of candidate regions may be refined by one or moreautomated decision-making steps. Both neural networks [93, 59] andrule-based [8] systems have been used in this capacity. Inputs to thesedecision-refinement stages come from one or more shape, edge, orhistogram features.

One of the most important issues regarding the use of CAD systems fornodule detection is the effective detection rate. Most computer-aidedschemes for the analysis of chest radiographs unfortunately haverelatively low sensitivity. In addition there is frequently a highfalse-positive rate (very low specificity) in these systems. A typicalresult is 2-7 false-positives per image, with sufficient sensitivity(over 70%) [93, 8]. Thus, the choice of operating point on the receiveroperating characteristic (ROC) curve for a particular CXR analysistechnique may be quite challenging.

In addition to the detection of pulmonary nodules, computer techniqueshave recently been studied for the classification of nodules found inchest radiographs. Computation of density gradients as a measure ofshape and surface irregularity has been used to distinguish benign frommalignant nodules [71]. Fractal texture analysis has also been appliedto the nodule classification problem with encouraging results [58, 87].

When assessing nodule growth for the prediction of malignancy, two ormore CT studies, separated by a suitable interscan interval, are needed.It is theorized that it may also be possible to estimate the probabilityof nodule malignancy based on size, shape, and density parametersdetermined from a single exam.

In a screening setting, 2D features are frequently measured to refinethe set of nodule candidates for detection. In this capacity, Giger etal. [22] used measures including perimeter, area, compactness,circularity, elongation, and location within the lung. Similarly,Toshioka et al. [84] used area, thickness, circularity, density mean andvariation, and location. For small nodule characterization, however,HRCT studies are required (see FIG. 3) and methods for high-resolutioncharacterization are being explored [42, 63, 51].

McNitt-Gray et al. [51] have explored two-dimensional (2D) shape metricsincluding aspect ratios, circularity, and compactness, as well as avariety of texture measures (e.g. correlation texture and differenceentropy). Their study also employed bending energy as a measure ofsurface characteristics. Three-dimensional measures of surfacecharacteristics were studied by Kawata et al. (42). Using dataisotropically resampled to 0.33 mm³ voxels, they measured Gaussian andmean curvature as metrics of irregularity in nodule margin.

Another challenging problem in the study of pulmonary nodules is theanalysis of non-solid nodules, or ground-glass opacities (GGOs) [18,11]. These lesions are difficult to characterize using traditional sizeand shape metrics.

Patient demographic information has been used as an adjunct toradiographic analysis in the prediction of small pulmonary nodulestatus. Several patient characteristics have been shown to be highlycorrelated with nodule malignancy. These include advanced age andsmoking history (pack-years) [81, 82, 47].

To study the pulmonary nodule using computer-aided techniques, we mustformulate one or more models of what a pulmonary nodule is, inparticular, with respect to its radiographic appearance. Pulmonarynodules are lung neoplasms that fall into two general categories: benignand malignant. The separation between these categories is generally onebased on rate of growth. Most benign lesions do not grow, or grow muchmore slowly than malignancies. There are, however, some notableexceptions. Infections in the lung may exhibit extremely fast growth.Also, it has been suggested that there may be some indolent,slow-growing malignant species. These malignancies are by far theexception, however, and doubling times for most malignant nodules rangefrom 30-400 days [47].

Pulmonary nodules appear radiographically as round, opaque lesions withdensity slightly more dense than that of water (˜0-100 HU). Theirposition in the lung and immediate surroundings, from an image-analyticperspective, differentiates them into one of the following categories:(a) well-circumscribed—the nodule is located centrally in the lung,without significant connections to vasculature; (b) vascularized—thenodule is located centrally in the lung, but has significantvascularization (connections to neighboring vessels); (c) pleuraltail—the nodule is near the pleural surface, connected by a thinstructure (“pleural tail”); and (d) juxtapleural—a significant amount ofthe nodule periphery is connected to the pleural surface.

The above categories describe the types of solid pulmonary nodules seenin CT images. A second class, known currently as ground glass opacities(GGO) have a significantly less uniform radiographic appearance and willbe considered separately later in this work.

One of the best predictors of nodule malignancy is growth rate.Similarly, nodule size, is also highly correlated with malignancy [82].Significantly larger nodules (>3 cm in diameter) are more likelyrepresent lung cancer [82]. Overall size has therefore been used incharacterization and classification studies [26, 82].

In addition to these benefits, volumetric measurement allows for theidentification of anisotropic growth. It is commonly assumed thatpulmonary nodules grow uniformly in all dimensions.

In 1956, Collins et al. [12] succinctly described the relationship ofgrowth and cancer, as well as the need for quantitative measurements oftumor growth:

-   “The definition of cancer, its diagnosis and its prognosis all    depend upon description of growth. To the layman a synonym for    cancer is a ‘growth.’ There are no quantitative terms for the    description of growth or growth rate in clinical use.”    Their seminal paper continues, providing some of the most important    early discussion of tumor “doubling time,” a measure of the amount    of time (usually expressed in days) for a tumor to double in volume.    Nodule growth, and the rate at which it takes place, is perhaps the    feature most predictive of malignancy. Nodule doubling time and    calcification in “benign” patterns have been the only universally    accepted diagnostic tools of thoracic radiologists in the    differentiation of benign from malignant lesions.

Traditional measures of the size and shape of pulmonary nodules weremade manually on chest radiographs using calipers. Nodule size wascharacterized by a measure of diameter (based on the assumption thatnodules appeared roughly circular), or as a combination of several“diameter” measurements.

Assessment of pulmonary nodule shape has been based on subjective shapecharacteristics including notions of roundness, spiculation, lobulation,and sharpness of edge. Nodule surface characteristics are frequentlyused to help classify benign from malignant nodules. In particular, thepresence of spiculation on the surface, as well as lobulation of theoverall nodule shape, have been reported as more common in malignantnodules [47, 82]. As with measurement of nodule size, there have beendifficulties in making reproducible shape characterizations, due to thelack of a precise mathematical basis and measurement technique.

Several techniques for measuring nodule surface characteristics havebeen described in the literature. Huo et al. [32, 33] used analysis ofradial edge gradients (a technique originally developed for mammography)in two-dimensional CT images. Kawata et al. [42] used a method ofcurvature estimation based on 3D image gradient analysis. Suchtechniques for curvature estimation, and the prerequisite 3D gradientanalysis and edge detection have been described in detail [5, 54].

Detection of pulmonary nodules in screening CT scans is challenging dueto the desire to limit radiation dose and scan time, both of whichresult in a reduction of information. The need to scan the full lungvolume in a single breath hold motivates the use of large intervalsbetween images (10 mm). An additional reduction in radiation doseachieved through low tube current further affects image quality. Withthese constraints, computer algorithms for the detection of pulmonarynodules have traditionally focused on two-dimensional (2D) algorithms orserial 2D methods that examine connectivity between nodule candidates inadjacent slices.

Automated CAD systems for the detection of pulmonary nodules have beenbased on a multistage approach, as illustrated in FIG. 8. The mainstages are: (a) preprocessing, (b) identification, and (c)classification. Preprocessing involves noise filtering and isolation ofthe lung region from the image data. Identification is the use of imagefilters and knowledge based reasoning to develop a set of nodulecandidates from the lung volume being searched. Classification is theanalysis of nodule candidates, normally as part of a rule-based decisionsystem, to reject other structures (e.g. vessels) from consideration andto assess the likelihood that each remaining candidate is a true nodule.The classification stage in detection systems is sometimes extended toinclude the separation of benign from malignant lesions, however, thismay only achieve limited results with the low-dose, low (axial)resolution studies common to initial screening exams. Characterizationand classification techniques that lead to prediction of malignancy aregenerally performed using high-resolution CT (HRCT).

A number of automated CAD systems for the detection of pulmonary nodulesin low-dose screening studies have been described in the literature [94,22, 84, 7, 4].

A typical lung cancer screening protocol consists of a low-dose (140kVp, 40 mA) CT scan of the full lung volume using 10 mm collimation,reconstructed at 5 mm intervals, with an in-plane resolution of 0.5 mmor better. Following reconstruction, the resultant voxels are highlyanisotropic (0.5×0.5×10 mm). Although the in-plane resolution may besufficient for the detection (and preliminary characterization) of manylesions, small pulmonary nodules (less than 10 mm in diameter) mayappear in only one to three contiguous images. This anisotropy presentsspecial challenges in distinguishing between small pulmonary nodules andfine vasculature.

Blood vessels may have either a linear or elliptical cross-section in asingle CT image, depending on orientation. The more perpendicular to theimage plane a vessel lies, the more circular the resultant image.Therefore, it is a non-trivial problem to distinguish between smallpulmonary nodules and small vascular cross-sections as the low axialresolution of screening studies limits the three-dimensional informationneeded to fully characterize these objects. The problem is furtherexacerbated by the fact that nodules are frequently attached to thepulmonary vasculature or to the pleural surface, increasing theirdifficulty of detection as they may be considered part of eitherconfounding structure due to the similar density values found in manynon-calcified nodules.

The preprocessing stage identifies the region of the CT images thatcorrespond to the lung volume. Algorithms have been developed toautomatically detect the thoracic cavity, lung volume, tracheobronchialtree, and pulmonary vasculature.

The first step in most pulmonary CAD systems is to identify the lungboundary, thereby eliminating other thoracic structures from the nodulesearch region. This is typically accomplished using serial 2D gray-levelthresholding techniques on each CT image. Spatial filtering may beapplied to reduce the effect of noise on the boundary segmentation. Oncea preliminary boundary has been determined, geometric filtering may alsobe used to eliminate discontinuities in the boundaries determined inadjacent slices, as well as to reduce the contribution of imageartifacts.

An example of this procedure is shown in FIGS. 4 through FIG. 6. FIG. 4shows a single 2D image from a low-dose screening scan. FIG. 5 shows thelung region in this cross-section produced using the techniquesdescribed above. A three-dimensional surface-shaded rendering can alsobe produced, as is shown in FIG. 6.

Giger et al. [22] used histogram analysis to determine appropriatethresholds for lung segmentation. Toshioka et al. [84] used splinerepresentation and curvature analysis to refine lung boundaries based onstandard thresholding as well as the detection of rib locations. Thisenables the detection of juxta-pleural nodules, which can be excludedfrom the lung volume using simpler threshold-based methods. Anatomicalmodels were developed by Brown et al. [7] as the basis for aknowledge-based system that utilized fuzzy set membership based ondensity, spatial locality, and position. These models were then used inthe determination of lung, tracheobronchial tree, and thoracic cavityregions. In our studies, we have found that simple thresholding combinewith geometric filtering has produced good results in identifying theseregions.

Another preprocessing step in pulmonary CAD systems is theidentification of the tracheobronchial tree and pulmonary vasculature.Segmentation of the airways may be accomplished by seeded region-growingtechniques. An example of this algorithm applied to screening data isshown in FIG. 12. It is evident that the low axial resolution ofanisotropic screening data limits the extent to which the bronchi can befollowed. Tree-based algorithms that explicitly track the bifurcationsof the bronchi and pulmonary vasculature have been described by Wood etal. [90]. Related methods based on morphological filtering [97] andfuzzy logic [103] as well as the benefit of airway subtraction to noduledetection have also been described [16].

Once preprocessing has eliminated regions outside the lung fromconsideration, the next step in nodule detection is the identificationof nodule candidates from the remaining lung volume. This is typicallydone through the use of one or more gray-level thresholding operations,as the nodules are of higher density than the surrounding lungparenchyma. Yamamoto et al. [94] combined single and double-thresholdingwith connected-component analysis to identify nodule candidates that arewithin a selected density range and size. In addition, morphologicalfilters, including disk and ring-shaped kernels as well as 2Dmorphological opening (“rolling ball filter”), were used to selectcandidates. Multiple thresholding based on histogram analysis has beenused to select nodule candidates and differentiate them from vessels[22]. In this system, while increasing the threshold, vasculature thatbifurcates along the scanner axis will produce one or more “daughternodes” in adjacent slices, thereby reducing the likelihood that thecandidate is a true nodule.

The final stage of a pulmonary nodule detection system, classification,is the refinement of a set of nodule candidates based on size, shape,and other measures. A likelihood estimate that the remaining candidatesare nodules is also often assigned. The techniques used in refining theset of nodule candidates are somewhat related to those used in thecharacterization and classification of nodules in high-resolutiondiagnostic studies. The difference, however, is that the aim of nodulecandidate classification is to separate nodules from non-nodules, whilecharacterization and classification in diagnostic HRCT is used toseparate benign from malignant nodules. With this distinction in mind,many of the nodule feature metrics used in detection and diagnostic CADsystems describe similar size and shape characteristics. It is importantto note, however, that the standard low-dose screening study issignificantly anisotropic, and therefore largely limits the feasibilityof performing measurements in three dimensions. For this reason,features used for nodule candidate classification have been primarilylimited to 2D metrics or analysis of connectivity between adjacent 2Dslices.

When a pulmonary nodule is detected at initial screening,high-resolution computed tomographic (HRCT) images may subsequently beacquired to further study the lesion. These studies typically involvescanning the region-of-interest (ROI) containing the nodule at a highin-plane resolution (0.3-0.5 mm), with small beam collimation (1.0 mm)and 1:1 pitch. At the time of HRCT, a diagnostic CT scan of the chest isalso commonly done for use in patient management. FIG. 3 shows acomparison between low-dose screening, diagnostic, and HRCT studies.

An important consideration for HRCT studies is that the entire 3D volumeof the lesion should be acquired without significant artifact (due torespiration or other patient motion). In traditional high-resolutionstudies, only the several images containing the largest 2D nodulecross-section were considered vital, as measurements would be madevisually or with calipers. Thus, it is an object of the presentinvention to provide the capability for enhanced size, shape, anddensitometric measurement of objects found in the body.

SUMMARY OF THE INVENTION

The present invention includes methods, systems, and computer readablemedium for producing and analyzing three-dimensional images. Inparticular, the present invention relates to image representation andanalysis for diagnosis and possible treatment of humans.

With respect to image production, the present invention includesproviding a three-dimensional image representation of an object whichdetectably attenuates a signal resulting from scanning. A signalcorresponding to the three-dimensional image is obtained and can besupersampled in three dimensions to a degree greater than the highestspatial frequency of the signal.

Scanning-medium used herein can include x-rays, magnetic field (asprovided by MRI), ultrasound, laser, electrical current, visible light,ultraviolet light, infrared light, and radio frequency. A preferredembodiment of the present invention contemplates the use of x-rays suchthat the signals produced therefrom provides CT tomography forprocessing and analysis.

Preferably, supersampling includes imposing three-dimensional units ofimage representation (voxels) on a three-dimensional image signal andincreasing the number of such units in order to reduce the contributionof error per unit. Preferably, the supersampling is conducted inthree-dimensional isotropic space.

This method can be conducted on a host body, especially a patient,wherein the object is a tissue mass such as a growth, e.g., a tumor. Ithas been found especially useful for detecting and analyzing pulmonarynodules, and, especially, small pulmonary nodules which are equal to orless than 10 mm in diameter.

After the image has been supersampled, it can be further processed bysegmenting the object from other structures detected in the signalretrieved from the scan. Such further processing can include subtractingattachments such as vascular attachments from the object, and,subtracting the object from adjacent structures such as pleuralstructures. The image can be further processed by restoring volumetriccharacteristics of the object.

Another aspect of the present invention includes three-dimensionalsegmentation of a region of interest and a host body to differentiate anobject found in the region. The region of interest is scanned with, forexample, x-ray, to provide a signal corresponding to a three-dimensionalrepresentation of the region upon which three-dimensional units ofrepresentation are imposed. The units are then identified as havingobject signal, background signal and as boundary units which have bothobject and background signal. The image representation units are thensubjected to thresholding in order to separate background from object,followed by three-dimensional connected component analysis to identifythose images as objects contiguous in three-dimensional space.

The connected component analysis can include component labeling in threedimensions such as by recursive connected component labeling oriterative connected component labeling.

Furthermore, the segmentation procedure can include the step ofthree-dimensional supersampling of the three-dimensional signal image.

Segmenting can also include morphologically filtering images identifiedas a result of the steps set forth above. Such filtering can includethree-dimensional opening of the images where by small connectionsbetween regions are broken and sharp edges are removed, such openingincluding erosion followed by dilation of the images in threedimensions. The method can also include a closing of thethree-dimensional region of the background signal resulting from theprevious steps, such closing including dilation followed by erosion ofthe background in three dimensions.

The segmentation can also include three-dimensionally regrowing imagesto approximate true surface characteristic of the images. In oneembodiment the regrowing is iterative constrained dilation in threedimensions.

The segmentation procedure can also include eliminating structuresadjacent images of interest, and, in the case of a nodule, suchstructures can include thoracic structures such as a pleural surface. Inorder to remove thoracic structures, the present invention contemplatesdetermining in three dimensions angles describing orientation of thesurface of a structure to be eliminated followed by performing anopening operation (based on the angles previously found) to detect amajority of the structure to the exclusion of the image of interest.Thereafter, the structure is three-dimensionally subtracted from theimage of interest. In a preferred embodiment, the step of determining inthree dimensions the angles describing orientation of the surface of thestructure can be conducted by three-dimensional moment analysis.Moreover, the method for eliminating adjacent structure can also includemorphologically filtering in three dimensions the image of interest.Finally, the signal can be processed to provide a smooth surfacerepresentation of the image of interest.

A smooth surface representation of a three-dimensional voxelrepresentation of a three-dimensional image can include: segmenting theimage to provide a segmented voxel three-dimensional signal; followed bymodified tessellating the three-dimensional image to provide athree-dimensional triangulated curve surface estimation; after which thetriangulated curve surface estimation can be smoothed by filtering.

Another aspect of the present invention is the ability to provide anestimated volumetric change time (DT) of an object which undergoeschange in size in the body, such as a growth or even an organ found in amammalian host, e.g., a human. This method contemplates obtaining ameasurement of a change of a volumetric characteristic of the objectover a period of time. The “change” referred to herein can be growth orregression of growth. The change of the volumetric characteristic isthen related to the time period of such change by an exponential changemodel whereby an estimated doubling time can be derived by comparing thevolume change over such time period. An “exponential change model” asused herein includes “exponential growth model.” Thus, when the term“exponential growth model” is used in the description, exponentialchange model is deemed to be implicit.

This method is implemented by obtaining two three-dimensional volumetricassessments of an existing object at a first time t₁ and at a latersecond time t₂. The volumetric assessment can be a volume V₁ followed bymeasurement of a volume V₂. Alternatively, the volumetric characteristiccan be a measure of the mean density at time t₁ and a measure of themean density at time t₂. Another embodiment contemplates combing achange retardation factor, R, with the standard exponential growthmodel.

In a preferred method of estimating the doubling time, the measurementscan be obtained from three-dimensional images provided from signalshaving three-dimensional data resulting from scanning the host, and, inother preferred embodiments, when such image has been supersampled,segmented, and both supersampled and segmented. This method isparticularly useful in conducting diagnosis and possible treatment ofpulmonary nodules, especially small pulmonary nodules which are usuallynot greater than about 10 mm in diameter.

An alternative embodiment of this method and system include the abilityto estimate doubling time of a an object which does not appear in afirst scan of the region of interest but does appear at a second scan attime t₂. In this case a minimum detectable change dimension β isidentified followed by calculating doubling time using a volumetricmeasurement at t₂ and a volume measurement derived from the minimaldetectable change dimension β over the time period between t₁ and t₂.

Another embodiment of the present invention related to time betweenscans includes a method and system for determining time between scanningexaminations, or interscan intervals. In this aspect of the invention, amaximum magnitude of error, ε, between volumetric measurements fromscanning and actual object size is developed. A minimumreliably-detectable percent volume change (α) is then identified basedon the maximum magnitude of error. A relationship is established betweenthe magnitude of error and object size, preferably a linearrelationship, and an expression using an exponential change model isderived in terms of change in volume and time. Based on the expressionresulting from the exponential change model relationship, a timeinterval to observe a percent change is determined. Preferably, thismethod and system is highly applicable to potentially cancerous growths,especially when the interscan interval is sufficiently short in durationto increase probability of early detection as cancerous. Once again thisprocedure is especially helpful in detection, analysis, and diagnosis ofsmall pulmonary nodules, e.g., those not greater than 10 mm in diameter.

The present invention also contemplates a method and system ofdetermining volumetric change in size over a period of time. In thismethod a doubling time of change is first determined, followed byestablishing a relationship based on a exponential change model of saidperiod of time between volumetric change time and a relative change involume, followed by determining the volumetric change over such periodof time.

Another method and system of analyzing three-dimensionalscanned-produced images includes three-dimensionally registering aregion of interest taken at separate times. In one method theregistration is effected by locating the center of mass (COM) in theregion of interest and aligning subsequent images of the region ofinterest for analysis. This can be accomplished by iterativelydetermining the center of mass as explained in detail hereinafter.

In another manifestation of this embodiment, the three-dimensionalregistration of regions of interest is conducted by three-dimensionalcorrelation. In this case two images of region of interest are selectedwherein the limits of search in such region are specified. The secondregion of interest is translated to all locations in the search regionand the three-dimension correlation between the two regions of interestare determined to find the greater value of a match metric.

Analysis of three-dimensional registration of two regions of interestcan also include a process of weighting such as by use of Gaussianwindow, a triangular radial window, or a rectangular radial window,which can be used to reduce unwanted confounding structures found in theregion of interest. The three-dimensional weighting of the region ofinterest can also include method of moment analysis to provide densitydistribution of the weighted region of interest.

Yet another aspect of the present invention is a method of estimatingcurvature of an object which has been scanned to provide a signal havinga three-dimensional representation of such object. This method includesproviding a three-dimensional triangularly tessellated representation ofthe object followed by determining the surface normal to all trianglesin the tessellated representation. In this regard it is noted that everyplane surface has a normal, thus each triangle has a surface normalwhich is determined in accordance with this method. Next, a surfacenormal at each vertex is then calculated based on the surface normalsresulting from the previous step. Then the angular differences betweenthe surface normal at each vertex (designated home vertex, V_(i)) andthe vertex normals of all vertices adjacent to the home vertex aredetermined. Finally, the curvature is estimated at each vertex of theobject based on the angular differences resulting from the previousstep. The relationships and algorithms supporting this method have beenset forth in detail in the Detailed Description of the invention.

The present invention also includes a method of providing a curvaturemetric analysis of a three-dimensional object which has been scanned toprovide a single corresponding to such object. In this method, ameasurable three-dimensional representation of an object is providedwhich includes a surface. Surface curvature is estimated of such surfaceat a number of discreet points on the surface to provide a number ofsurface curvature estimates. Then a frequency distribution of thesurface curvature estimates resulting from the previous step isdetermined. Subsequent to determining the frequency distribution of thesurface curvature estimates, various aspects of the distribution can bedetermined, including, but not limited to variance, range ofdistribution, mean of the distribution, minimum of the distribution, andmaximum of the distribution. Preferably the method set forth above canbe applied to an object which is a growth in a human such as a pulmonarynodule, especially one which is not greater than 10 mm in diameter.

Moreover, the analysis based on three-dimensional registration can alsoinclude estimated doubling time of an object based on a change in themean density of the region of interest. In this case, an object such asan object image, is Gaussian-weighted and the estimation can be aniterative technique to select an appropriate standard deviation, σ, forweighting. The doubling time has been calculated based on standardizedmean density and a stopping criteria, ε.

As a result of the present invention, three-dimensional resamplingtechniques have been developed and evaluated which mitigate partialvolume effects in the assessment of object volume. There has been asignificant improvement in the reproducibility of volumetric estimatesafforded by the techniques of the present invention over use of standardreconstructed CT data. For example, the double time estimates areconsistent with pathologic diagnosis and are more appropriate estimatesthan those based on 2D metrics. Furthermore, as a result of estimatingcurvature based on a smoothed tessellated polygonal surface model, theactual nodule surface is approximated given the available resolution.

Another benefit of the present invention is to establish new techniquesfor computer-aided diagnosis (CAD) of small pulmonary nodules (SPNs). Acollection of three-dimensional size, shape, and density metrics havebeen developed for characterization of small pulmonary nodules. Thisdevelopment has enhanced the ability to differentiate benign frommalignant lesions.

Furthermore, new methods have been developed for possible computer-aidedanalysis and diagnosis of small pulmonary nodules in computedtomographic data. Small pulmonary nodules, those with a diameter of lessthan 10 mm in diameter, are now commonly detected and analyzed accordingto their size, shape, density and growth characteristics.

For a better understanding of the present invention, together with otherand further objects, references made to the following description takenin conjunction with the accompanying drawings and the scope of theinvention will be pointed out in claims.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is an illustration of Beer's Law in a homogeneous material;

FIG. 2 is an illustration of Beer's Law in a heterogeneous material;

FIG. 3 illustrates a comparison of screening and diagnostic CTprotocols;

FIG. 4 illustrates a single slice from a whole-lung CT scan;

FIG. 5 illustrates the lung region identified from the CT slice shown inFIG. 4;

FIG. 6 illustrates a 3D reconstruction of the lung surface;

FIG. 7 illustrates a 3D reconstruction of the trachea and main bronchi;

FIG. 8 illustrates a multi-stage approach for the automated detection ofpulmonary nodules in screening scans;

FIG. 9 illustrates a well-circumscribed small pulmunary nodule;

FIG. 10 illustrates a small pulmonary nodule with confoundingvasculature;

FIG. 11 illustrates a small pulmonary nodule with a pleural tail.

FIG. 12 illustrates a small juxtapleural nodule;

FIG. 13 illustrates a model of a vascularized nodule;

FIG. 14 illustrates a volume overestimation at point of vesselattachment;

FIG. 15 illustrates a volume of a spherical cap;

FIG. 16 illustrates a volume overestimation of an 8 mm spherical nodulewith vascular attachments of varying sizes following opening with aspherical kernel;

FIG. 17 illustrates a model of a juxtapleural nodule;

FIG. 18 illustrates a two-dimensional sampling of a circle where nequals 0;

FIG. 19 illustrates a two-dimensional sampling of a segment of a circlewhere n equals 1 and 2;

FIG. 20 illustrates a two-dimensional area estimation of a segment of acircle;

FIG. 21 illustrates a two-dimensional image of a small pulmonary nodule;

FIG. 22 illustrates a histogram of gray-level intensities in noduleimage;

FIG. 23 illustrates a two-dimensional segmentation of a small pulmonarynodule using iterative threshold selection;

FIG. 24 illustrates segmentation results by varying threshold from −800HU to 100 HU;

FIG. 25 illustrates examples of morphological operations using aspherical kernel;

FIG. 26 illustrates a 2D morphological opening on a pulmonary nodule;

FIG. 27 illustrates a three-dimensional morphological opening forvascular subtraction;

FIG. 28 illustrates a three-dimensional iterative morphologicalfiltering for vascular subtraction;

FIG. 29 illustrates a three-dimensional iterative morphologicalfiltering for vascular subtraction.

FIG. 30 illustrates segmentation of a synthetic vascularized nodule;

FIG. 31 illustrates sensitivity of vascular subtraction as a function ofvessel diameter and kernel size;

FIG. 32 illustrates sensitivity of vascular subtraction to kernel sizein segmenting in vivo pulmonary nodules;

FIG. 33 illustrates an example of the pleural surface removal technique;

FIG. 34 illustrates an example of the pleural surface removal technique;

FIG. 35 illustrates a 15 conical 2×2×2 neighborhoods and associatedtriangulations in the marching cubes algorithm;

FIG. 36 illustrates a comparison of 3D surface representations of aknown sphere;

FIG. 37 illustrates a comparison of exponential nodule growth for threedoubling times TD=75 (aggressively malignant), TD=150 (malignant), andTD=500 (benign);

FIG. 38 illustrates RMS percent error in 35 volume measurements ofspherical silicone nodules;

FIG. 39 illustrates CT images of a small nodule at baseline and 33 dayslater;

FIG. 40 illustrates 3D visualization of nodule at baseline and 33 dayslater;

FIG. 41 illustrates two-dimensional images and segmentation of amalignant small pulmonary nodule at baseline and 36 days later;

FIG. 42 illustrates 3D visualization of malignant nodule exhibitingnon-uniform growth;

FIG. 43 illustrates an ellipsoid with its principal axes labeled;

FIG. 44 illustrates a silhouette of an object whose shape will bediscretized by the sampling procedure;

FIG. 45 illustrates the discretized representation of the object shownin FIG. 44;

FIG. 46 illustrates a polyline representation of the object shown inFIG. 44;

FIG. 47 illustrates a filtered (smooth) polyline representation of theobject shown in FIG. 44;

FIG. 48 illustrates a comparison of 3D surface representations of asmall pulmonary nodule;

FIG. 49 is an illustrative patch of a 3D tessellated surface;

FIG. 50 illustrates a computation of the surface normal of eachtriangle;

FIG. 51 illustrates a calculation of the surface normal at each vertex;

FIG. 52 illustrates surface normals forming the basis of curvatureestimation at a vertex;

FIG. 53 illustrates an atomic measure of curvature based on the angulardifference between vertex surface normals;

FIG. 54 illustrates a reduction in variance of curvature estimates ofsynthetic spheres when using average normal curvature over meancurvature;

FIG. 55 illustrates a surface curvature analysis of a malignant smallpulmonary nodule;

FIG. 56 illustrates a mean of average normal curvature distribution insynthetic spheres;

FIG. 57 illustrates a spherical perturbation model for surface curvatureanalysis;

FIG. 58 illustrates a frequency distribution of curvature estimates in asynthetic sphere;

FIG. 59 illustrates a frequencey distribution of curvature estimates insurface perturbations;

FIG. 60 illustrates a frequency distribution of curvature estimates of aspherically-perturbed synthetic sphere;

FIG. 61 illustrates a frequency distribution of curvature estimates of aspherically-perturbed synthetic sphere;

FIG. 62 illustrates a variance of curvature estimates as a function ofsmoothing;

FIG. 63 illustrates a Receiver-Operating Characteristic (ROC) curve forthe complete set;

FIG. 64 illustrates a Receiver-Operating Characteristic (ROC) curveusing cross-validation;

FIG. 65 illustrates an iterative center of mass determination in a CTscan of a small pulmonary nodule; and

FIG. 66 illustrates a correlation-based registration of nodule ROIs.

DETAILED DESCRIPTION OF THE INVENTION

The present invention relates to methods and systems for conductingthree-dimensional image analysis and diagnosis and possible treatmentmodalities relating thereto. For the purpose of disclosing the presentinvention and the present preferred embodiments thereof, the detaileddescription has been broken down by subject matter in order to moreclearly explicate the various aspects of the invention. Headings havebeen provided for the sake of convenience and focus of subject matter,but are not meant in any way to limit the scope of applicability of thesubject matter found in one section. Indeed, each technical aspect isgenerally considered germane to subject matter found in other sections.Thus, one must read all of the sections as related to the extent thatthe technology disclosed therein can be applied to the technology foundin other sections.

Furthermore in this regard, a certain degree of redundancy has beenbuilt into the detailed description in order to provide as much aspossible overview and correlation of the technology between sectionsfound herein.

A system in accordance with one embodiment of the present invention caninclude a scanner, a processor, a memory, a display device, inputdevices, such as a mouse and keyboard, and an internal bus connectingthe various components together. The system can be coupled to acommunication medium, such as a modem connected to a phone line or awireless network, or an Internet connection. Since the components of acomputer and how they operate are well known to one of ordinary skill inthe art, they will not be discussed here.

A scanner includes any device capable of imposing a scanning medium onsubject, e.g., a patient, to produce a signal attenuated by objects inthe field of the scan. These devices include, but are not limited to, CTscanners, MRI, PET, ultrasound, optical tomography, and electricalimpedance tomography. A preferred scanner is a helical CT scanner whichwas used in the work and analysis reported herein.

A processor as referred to herein is a system of one or more electroniccomponents which are statically or dynamically configured (in hardwareor software) to implement an algorithm or process. Examples includesystems based on microprocessors, application-specific integratedcircuits (ASICs), field-programmable gate arrays (FPGAs), and/ordiscrete components. Any type of storage device, such as a RAM or ROM,can be used as the memory in the system. Additionally, the memory can bea magnetic, optical, magneto-optical, or solid state drive that iscoupled to the processor and which can receive portable memory devices,such as floppy disks, hard disks, optical disks, or solid-state memoryon which the method of processing computed tomography images inaccordance with the present invention is programmed. All media which canstore information, operating instructions, etc., and can be read by acomputer for processing, storage, etc. is referred to collectivelyherein as “machine readable-medium.” The method of processing computedtomography images in accordance with the present invention can be storedas an executable program in the memory.

Section 1: Computer Assisted Diagnosis (CAD)—Virtually all aspects ofthe present invention can be used to provide computer assistance in theprocess of diagnosis for detecting, treating, and monitoring thebiological condition(s) of a patient. The degree to which the techniquesand technology disclosed herein are used and/or relied on depends on thecondition, the patient, the predictability of outcome, the physician andother factors too numerous to detail herein. To the extent thetechnology disclosed herein can be considered as relevant to (e.g., as abasis for or as part of) Computer Assisted Diagnosis (CAD), each aspectof the invention claimed herein includes such relevancy to CADapparatus, systems and methods.

Section 2: Nodule Size and Shape—As a means for explanation andanalysis, a basic model for the appearance of a small solid pulmonarynodule is defined as roughly spherical, with a diameter less than 1 cm.The nodule has a density (CT attenuation) significantly higher than thatof the surrounding lung parenchyma. The nodule form may be confounded byother neighboring/attached structures, including vessels and the pleuralsurface. As the concern is the analytic characterization of the noduleitself, an important task is image segmentation the delineation of thenodule boundary and removal of other structures from consideration. Theimage segmentation process is one of the most challenging aspects of anycomputer vision task. Segmentation methods for pulmonary nodules in CTscans are discussed herein. Segmentation can also be viewed aspreprocessing of the image region of interest (ROI) to yield athree-dimensional representation of the nodule which conforms to ourbasic model of a pulmonary nodule.

Consider this basic model in the case where the nodule is wellcircumscribed. An image of a well-circumscribed pulmonary nodule isshown in FIG. 9. In this case, the segmentation may be accomplishedusing a simple gray-level thresholding operation. Selection of theappropriate threshold involves the following issues. First, a thresholdcan be selected that best separates the nodule (foreground) from thelung parenchyma (background), which in the well-circumscribed case wouldhave a bimodal intensity histogram. This separation may be achieved byselection of a threshold that falls at the midpoint between the peaks ofeach mode of the histogram, or by using a process (and relatedalgorithmic relationship) that considers the area under the histogram ineach region. Thresholding methods are discussed in Section 4.

The second consideration when selecting a threshold is in evaluating itschange in size over time. Thus a single threshold should be selectedwhich is appropriate across multiple studies of the same nodule. Inaddition, a standard threshold may be determined for a general class ofwell-circumscribed nodules such that the size of nodules can be comparedacross time and multiple subjects.

Confounding Structures—Well-circumscribed pulmonary nodules with theircharacteristically high contrast to the surrounding lung parenchyma areunfortunately not the only class of pulmonary nodules commonly seen inCT studies. Measurement of the size and shape of pulmonary nodules isfrequently confounded by the presence of other structures nearby oradjoining the nodule. A simple two-level model based on CT attenuationis insufficient to completely separate the pulmonary nodule underconsideration from these other structures, as they frequently exhibitsimilar density characteristics.

Vasculature—Pulmonary nodules, as other cellular structures in the body,require oxygen and nutrients to sustain growth. These are supplied tothe nodule by the local vasculature (blood vessels). The radiographicappearance of these vessels ranges from undetectable in peripherallesions, to significant in more central ones. FIG. 10 shows an image ofa small pulmonary nodule with surrounding and attached vessels.

A variety of physiological factors underly the presence of vascularattachments to pulmonary nodules. In the context of this work, however,it is sufficient to state that vessels may appear nearby orsignificantly attached to nodules, which complicates the nodulesegmentation problem. Method (and related algorithms) for thesubtraction of blood vessels during image segmentation are discussed inSection 4.

Pleural Surface—Peripheral pulmonary nodules often exhibit some degreeof attachment to the pleural surface (on the periphery of the lung,compressed against the external boundary of the thorax). One possibletype of attachment is the so-called “pleural tail sign.” This appears asa small, thin structure of nodule-like attenuation connecting the noduleto the pleural surface. A nodule with a pleural tail is shown in FIG.11. In such cases, segmentation is relatively straightforward andmethods similar to those used to remove confounding vessels may beemployed.

A second class of pleurally-attached nodules are those that arejuxtapleural. These nodules share a significant amount of their surfacewith the pleura. For this reason, delineation of the boundary betweenpleura and nodule can be quite challenging. A small juxtapleural noduleis shown in FIG. 12. Segmentation techniques for juxtapleural nodulesare described in Section 4.

Nodule Models—Segmentation and measurement of pulmonary nodules in CTdata benefits from a more formal mathematical model and description ofeach nodule type we wish to describe. The four basic solid nodule typesare as follows: well-circumscribed, vascularized, pleural-tail, andjuxtapleural. Having given a qualitative description of each one ofthese types, we may now formalize mathematical models for each.

Well-Circumscribed—A well-circumscribed nodule is one that is locatedcentrally in the lung, without significant connections to vasculature. Asimple model for the nodule would be based on the fact that there iscontrast in x-ray attenuation and, therefore, voxel intensity, betweenthe nodule and the surrounding lung parenchyma. This model a set N ofvoxels in the nodule would then include all voxels greater than or equalto a given intensity threshold:N={v(x,y,z)|v(x,y,z)≧T}  (2.1)

This model would be very sensitive to noise, however, as any voxel inthe region of interest meeting the threshold criterion would beincluded. A better model is to include connectivity between the voxelsin the nodule. For each voxel in the set, then, there would have to beat least one neighboring voxel also in the set. Note that weintentionally exclude the trivial case where the nodule would be asingle voxel. This connected nodule model could be defined in terms ofthe the set of voxels adjacent to a given voxel. This set will bedenoted as adj(v(x, y, z)). Voxel connectivity will be described belowin detail. Our model for the well-circumscribed nodule is then the setof connected voxels greater than or equal to a given intensitythreshold:N _(c) ={v(x,y,z)|(v(x,y,z)≧T)

(∃n|nεadj(v(x,y,z)),nεN _(c))}  (2.2)

There is the possibility that there may be more than one connectedcomponent A in the region of interest (ROI). Considering the fact thatthere may be more than one set of connected voxels meeting the thresholdcriterion, we will elect to choose the largest of these sets, N_(max),(that which has the largest volume):N_(max),=max A _(i)  (2.3)where each of the eligible sets is defined asA _(i) ={v(x,y,z)|(v(x,y,z)≧T)

(∃n|nεadj(v(x,y,z)),nεA _(i))}  (2.4)

Thus, we will assume that the nodule will be the largest connectedcomponent (meeting the threshold criterion) in the ROI.

Vascularized—A vascularized nodule is one that is located centrally inthe lung, but has significant vascularization (connections toneighboring vessels). Clearly, the use of the well-circumscribed nodulemodel for the segmentation and analysis of vascularized nodules isinappropriate, as that model of connectivity would include the vascularcomponents as part of the nodule. Instead, we define a model whichdistinguishes nodule volume from connected vessels geometrically.

Consider a model in which the vascularized nodule is the union ofspheres of diameter μ.N _(v) =∪s(x,y,z)|d _(s)=μ  (2.5)

Each of the voxels in each one of these spheres must satisfy a thresholdcriterion, as was true for the well-circumscribed nodule model, as wellas have no more than μ/2 distance from the sphere's center. Thefollowing equation describes such a sphere, s(x, y, z), centered atvoxel (x, y, z), where D is the Euclidean distance function between twovoxel locations.s(x,y,z)=∪{(v(x _(i) ,y _(i) ,z _(i))|(v(x _(i) ,y _(i) ,z _(i))≧T)

(D({x _(i) ,y _(i) ,z _(i) },{x,y,z})≦μ/2)}  (2.6)

We will additionally define a maximum diameter, λ, of any vessels whichmay be connected to the nodule volume. Then, if we choose μ such thatμ>λ, we ensure that the majority of each vessel volume will not beincluded in the nodule volume.

This model of a vascularized nodule is illustrated in two dimensions inFIG. 13. The illustration on the left shows all of the nodule and vesselarea (volume in 3D) which meets the thresholding criterion. Theillustration on the right shows how the nodule volume can be describedas the union of circles (spheres in 3D) of diameter μ, while effectivelyignoring the majority of the vascular component. The model of avascularized nodule as a union of translated spheres is closely alliedwith the notion of the opening operation in mathematical morphology.This topic and it's practical relation to the vascularized nodule modelwill be discussed below in greater detail. One of the issues of interestis the overestimation of nodule volume at the point of vesselattachment, as illustrated in FIG. 14. Note how the spherical“structuring element” allows the inclusion of a small amount of thevessel volume (outside the spherical nodule volume). For an idealizedmodel of a spherical nodule with a vascular attachment of diameter λ, itis possible to formally state the extent of this volume overestimation.The volume in the extraneous region, V_(ext), is the difference betweenthe portion (spherical cap), {circumflex over (V)}μ, of the structuringelement that extends beyond the line of length λ, and the spherical cap,{circumflex over (V)}_(s), of the spherical nodule volume beyond thesame line. These two regions are illustrated (in 2D) in the right-handpart of FIG. 14. The extraneous region, V_(ext) is that between thedotted line (the boundary of the large sphere) and the externalboundary. The largest perpendicular distance in this cup-shaped volumeis shown as v. The volume {circumflex over (V)}μ is the region betweenthe straight line of length A and the external boundary, while thevolume {circumflex over (V)}_(s) is that between the straight line andthe dotted boundary of the large sphere. Therefore, the extraneous,cup-shaped, volume is simply their difference:V _(ext) ={circumflex over (V)}μ−{circumflex over (V)} _(s)  (2.7)

FIG. 15 illustrates the volume of a spherical cap. The volume as afunction of r, the radius of the sphere, and a, the radius of theintersecting disk is

$\begin{matrix}{{\hat{V}( {r,a} )} = {\frac{\pi}{3}( {r - \sqrt{r^{2} - a^{2}}} )^{2}( {{2r} - \sqrt{r^{2} - a^{2}}} )}} & (2.8)\end{matrix}$

Using Equation 2.8, we may substitute the appropriate values of r and ato derive the volumes of {circumflex over (V)}μ, and {circumflex over(V)}_(s), and subsequently, V_(ext). For the structuring sphere, ofradius μ/2, the volume of the spherical cap is

$\begin{matrix}{V_{\mu} = {\frac{\pi}{3}( {{\mu/2} - \sqrt{ {( {\mu/2} )^{2} - {\lambda/2}} )^{2}}} )^{2}( {\mu - \sqrt{ {\mu/2} )^{2} - ( {\lambda/2} )^{2}}} )}} & (2.9)\end{matrix}$and for the large sphere of radius r_(s), the volume of thecorresponding spherical cap is

$\begin{matrix}{{\hat{V}}_{s} = {\frac{\pi}{3}( {r_{s} - \sqrt{ {r_{s}^{2} - {\lambda/2}} )^{2}}} )^{2}( {{2r_{s}} - \sqrt{r_{s}^{2} - ( {\lambda/2} )^{2}}} )}} & (2.10)\end{matrix}$The volume of the extraneous cup-shaped region, therefore, can bedetermined using Equations 2.7, 2.9, and 2.10. Note that theseexpressions are general for any sized spherical nodules, structuringspheres, and vessel diameters. More specifically, for a particularspherical nodule size and structuring sphere, the contribution ofmultiple vessels of different diameters can be determined. An example isillustrated in FIG. 16, where a spherical nodule of 8 mm in diameter(r=4) has a vascular attachment of varying diameter (λ). Each of thecurves in the graph corresponds to a λ value of 0.75, 1.0, or 1.25 mm.The graph illustrates the percent overestimation in nodule volume whenthe diameter of the structuring kernel used (μ.) is varied from 1.25001to 2.5 mm. In the worst case (λ=1.25 mm, μ=1.25001 mm), where thestructuring kernel is barely large enough to detach the vessel, thevolume overestimation is less than 0.19%. As the kernel diameter isincreased, this value decreases logarithmically.

Pleural Tail—The third variety of nodule is near the pleural surface andconnected to it by a thin structure, called a “pleural tail.” While itmay be advantageous in some respects to treat this group differentlythan those in all other groups, it may be modeled as one of the othertwo groups, as a vascularized nodule, or as a juxtapleural nodule. Inthe former case, consider a model similar to that of the vascularizednodule, where the nodule region is the union of spheres of diameter μ,wherein each voxel meets a fundamental threshold criterionN _(p) =∪s(x,y,z)|d _(s)=μ  (2.11)This expression is sufficient to exclude all of the pleural tail exceptfor a small, cup-shaped region. A full discussion was given in thecontext of vascularized nodules set forth above. When the volume of thepleural tail is to be considered jointly with the rest of the nodule,the juxtapleural nodule model tray be used, described in the nextsection.

Juxtapleural—In juxtapleural nodules, a significant amount of the noduleperiphery is connected to the pleural surface. FIG. 17 illustrates thebasic model. Four regions are defined: N, the nodule region; LP, thesurrounding lung parenchyma; PS, the pleural surface; and TW, thethoracic wall region. In CT images, the nodule region and the thoracicwall region both exhibit high-attenuation characteristics. Therefore,while simple thresholding is able to separate the nodule from the lungparenchyma, it is insufficient to separate the nodule from the outerthoracic region. In FIG. 17, a dashed line between points a and b hasbeen included to indicate one possible boundary between nodule andthoracic wall regions. Physiologically, the pleural surface separatesthe lung volume (containing the parenchyma) and the thoracic wall. Whena juxtapleural nodule is present, however, there may or may not beinvasion of the pleura. Whether the nodule has invaded the pleura or ismerely adjacent to it, voxel density (intensity) information is rarelysufficient to distinguish these regions.

In our model, we would like to consider the nodule region to be onlythat volume which is unlikely to be part of the thoracic wall. This canbe defined geometrically. We may define points a and b in thetwo-dimensional cross-section shown in FIG. 17 to be those points ofmaximum concavity for a given scale. The notion of scale in theestimation of concavity here is quite important, as the pleural surfaceboundary is likely to have local concavities due to image samplinggeometry and also physiologic/anatomic variation. Techniques fordetermination of the boundary between the nodule region and thoracicwall will be discussed below.

Section 3: Data Acquisition—Data acquisition and management arenon-trivial aspects to lung cancer research using computed tomography.All CT data used in this work were collected as part of the IRB-approvedCornell Early Lung Cancer Action Project (ELCAP). The data were acquiredusing GE High-Speed and LightSpeed CT scanners using the protocolsdescribed the Background and shown in FIG. 3.

Each two-dimensional CT image is a 512×512 image with 12-bit pixels,providing a dynamic range of 4096 gray-levels. The pixels are backedinto two-byte words (16 bits), simplifying the data alignment forstorage and processing. A given CT examination may generate as many as200 of these images, given several series and associatedreconstructions. A quick calculation illustrates the enormity (bycurrent technological standards) of the storage problem.Size_(image)=512·512·2 bytes=512 KB  (3.1)Size_(exam)=200·512 KB=100 MB  (3.2)

While the raw data for most studies is on the order of 50 MB, additionalstorage must be available for processed images, results of imageanalyses, and special visualizations. Given the large data storage andmanagement requirements of this project, standardized mechanisms weredeveloped for image accrual, processing, and archiving.

DICOM—The Digital Imaging and Communications in Medicine (DICOM) 3.0[13] communication protocol and image storage standards can be used asthe interface to the system. These are well-defined standards whichallow the acquisition of data from a variety of sources, providing asystem and vendor independent interface to the analysis system.

VisionX—The VisionX system, developed at Cornell University, is acollection of computer vision and image processing libraries which havebeen used in a wide range of image analysis research. The VisionXdevelopment system was the basis for the method, system, and softwaredeveloped in this invention.

Software for the transmission of DICOM images to and from PACS (PictureArchiving and Communication Systems) systems and conversion to andVisionX format was developed. In VisionX format, the DICOM headerinformation, selectively blinded with respect to patient-identifyinginformation, is retained to facilitate retransmission of data and/oranalyses to standard PACS systems.

Isotropic Resampling—Most CT data are sampled in an anisotropic space,where the resolution in each of the three dimensions is not equal. Thein-plane (x-y) resolution on modern scanners is typically better than1.0 mm in both dimensions. The axial (z) resolution, however, may beanywhere from near 1 mm to as high as 20 mm. A typical high resolutionscan protocol provides 0.5 mm in-plane resolution with 1 mm axialreconstructions, producing voxels that are (0.5×0.5×1.0 mm). This posesseveral problems. First, although high-resolution 2D measurements can bemade in any one of the highly-resolved CT slices, 3D measurements arehindered by the less-resolved axial dimension. Special anisotropic imagefilters need to take into account the different axial resolution, whichbecomes computationally cumbersome. A second, more important problem toconsider is that of partial volume effects.

Partial Volume Effects—The partial volume problem describes thequantized, discrete spatial nature of the CT image. Consider thefollowing two-dimensional example. A unit circle (radius=1) is perfectlysampled on a regular grid. Pixels are set if a sufficient percentage oftheir area corresponds to the interior of the circle. We may assume thisthreshold to be 50% of a pixel's area. The resolution of the system, atthe coarsest setting is equal to the radius of the circle. This isillustrated graphically in FIG. 18. For simplicity, most of thefollowing discussion will concern only one segment of the circle, theupper left 90° segment, and the corresponding quadrant of the square inwhich the circle is inscribed. The quadrant may be divided at any numberof regular intervals to define an appropriate sampling grid, but forthis example, we will choose this division factor to be powers of twofor the number of divisions in the x and y dimensions, with n theappropriate exponent of 2. In FIG. 18, for example n is equal to zero,and the quadrant is divided into one large pixel, with area equal to 1.FIG. 19 illustrates the cases where n is equal to 1 and 2.

With this model in place, we may discuss the accuracy of areameasurement of the circular segment as a function of sampling interval(image resolution). The true In the area of the circular segment is

$\frac{\pi\; r^{2}}{4} = {{\pi/4} \approx 0.7854}$In the case where n=1, the area is estimated to be 1.0, as more than 50%of the pixel is covered by the circle. When n is increased to 1 and 2,the area estimates improve. In the first case, three of four pixels areset, yielding an area estimate of 0.75. For n=2, 13 of the 16 pixels aresufficiently filled by the circle, giving an estimate of 0.8125.

FIG. 20 illustrates the cases in which n is equal to 1, 2, 3, and 4,showing only the quadrant of interest. In each case, the pixels thatwould be considered members of the circle are shaded. It is obvious fromthe figures that as the degree of subdivision of the sampling gridincreases, the error in area estimation of the circle decreases. It isalso possible to produce bounds for this error as a function of thedegree of subdivision. Consider, however, the case where the unifiedquadrant is the original pixel size in the acquired image. We may thencast our discussion in terms of the degree of resampling of the imagespace. We may define the supersampling ratio, s, to be the number ofdivisions of the original pixel size in each of the x and y dimensions.The n-th term of the series describing s, then, is simplys_(n)=2^(n)  (3.3)Similarly, the number of pixels, p, in the quadrant is given byp_(n)=2^(2n)  (3.4)

Now, consider that the error in area estimation is the sum of errorsmade in representing the circle boundary, and that we therefore onlyneed to consider the pixels through which this boundary passes toevaluate the error. Each of these pixels will each contribute to eitheran overestimate or underestimate of the true area. For example, in theupper left illustration (n=1) in FIG. 20, the upper right and lower leftpixels are set, each overestimating the true area of the segment. Theupper left pixel however, is unset, contributing to a localunderestimation of the area. Again, it should be clear that only thosepixels through which the true boundary passes contribute to this error.Pixels that are definitively in the image foreground (set) or background(unset) will be correct. In the upper right illustration (n=2) in FIG.20, there are 8 definitive foreground pixels and one definitivebackground pixel. The remaining seven, those on the boundary, eachcontribute some positive or negative error to the overall area estimate.

It is possible to formally define the number of such boundary pixels inthis model, those for which the relative proportion of foreground willdetermine whether they are set. The number of these boundary pixels, b,in this quadrant is defined asb _(n)=2^(n+1)−1  (3.5)The first few terms in the series, {1, 3, 7, 15, 31}, corresponding to0≦n≦4, describe the numbers of upper left quadrant boundary pixels shownin FIGS. 18 and 20. For example, in the lower left illustration in FIG.20 (n=3), the supersampling ratio, s, is 8, the number of pixels, p, is64, and the number of boundary pixels, b, is 15.

The magnitude of the contribution each boundary pixel to the overallerror is bounded by 0.5 times the area of the pixel, as we chose 50%coverage as the threshold for determining whether a pixel is set.Therefore, a conservative upper bound on the error for estimating thearea of the circular segment, ε can be expressed as

$\begin{matrix}{\varepsilon = {b \cdot \frac{a}{2}}} & (3.6)\end{matrix}$where a is the area of a pixel. Given that this is a unit circle, andtherefore the quadrant has unit area, the area of each pixel is

$\begin{matrix}{a_{n} = {\frac{1}{p_{n}} = {\frac{1}{2^{2n}} = ( \frac{1}{2} )^{2n}}}} & (3.7)\end{matrix}$Thus, the error bound for a given value of n may be restated as theseries where

$\begin{matrix}{\varepsilon_{n} = {{b_{n} \cdot \frac{a_{n}}{2}} = {( {2^{n + 1} - 1} )( 2^{- {({{2n} + 1})}} )}}} & (3.8)\end{matrix}$and, after some manipulation,

$\begin{matrix}{\varepsilon_{n} = ( {\frac{1}{2^{n}}\frac{1}{- 2^{{2n} + 1}}} )} & (3.9)\end{matrix}$Each of the two terms in this series converges to zero, with the latterconverging much faster. Essentially, as n increases by one, the numberof boundary pixels basically doubles, while their size, and thus theirerror contribution, is divided by four. This leads to an error boundthat is approximately halved with each increase in n.With this two-dimensional example, we can see how perfect interpolation,or resampling, of the image space can reduce errors in size estimation.In three dimensions, isotropic resampling (supersampling to an isotropicspace) allows segmentation decisions to be made on a super-resolved, orsupersampled, grid, allowing more accurate, consistent boundarydecisions to be made. The intensity in each image voxel is interpolatedto estimate the intensities in each of the sub-voxels in thesupersampled space. Thus, each original voxel intensity value isresponsible for several in the new image space, mitigating partialvolume effects in the original data, as a more precise spatial locationof the desired gray-level transition (boundary) can be identified.

Supersampling—Three-dimensional supersampling is the process by whichthe image is given a new representation in units bearing a higherspatial frequency than those in the original image, resulting in animprovement in image resolution. The original image may be anisotropicor isotropic and supersampled to either a more refined anisotropic orisotropic space. The process involves computing new voxel values atevery location in the new, more refined coordinate space. Each of thesevalues may be computed using one of a variety of interpolation schemesincluding trilinear, tricubic, and B-spline interpolation.

Trilinear interpolation—In the present invention, supersampling wasconducted using trilinear interpolation of each voxel value at location(x_(n), y_(n), z_(n)), where each new voxel value can be defined interms of the voxels' eight nearest neighbors in the original space,v(x_(i), y_(j), z_(k)), where

$\quad\begin{matrix}\begin{matrix}{i = \lbrack 0 \middle| 1 \rbrack} \\{j = \lbrack 0 \middle| 1 \rbrack} \\{k = \lbrack 0 \middle| 1 \rbrack}\end{matrix} & (3.10)\end{matrix}$For notational simplicity, we define three distances between v(x_(n),y_(n), z_(n)) and v(x_(i), y_(j), z_(k)) in terms of the perpendiculardistance along each of the x, y, and z Cartesian axes.

$\begin{matrix}\begin{matrix}{{dx}_{i} = {{x_{n} - x_{i}}}} \\{{dy}_{i} = {{y_{n} - y_{j}}}} \\{{d\; z_{k}} = {{z_{n} - x_{k}}}}\end{matrix} & (3.11)\end{matrix}$We also define the resolution in each dimension in terms of the distancebetween adjacent voxels in each dimension.

$\begin{matrix}\begin{matrix}{{rx} = {{x_{1} - x_{0}}}} \\{{ry} = {{y_{1} - y_{0}}}} \\{{rz} = {{z_{1} - x_{0}}}}\end{matrix} & (3.12)\end{matrix}$Given the above, the value of the new voxel based on trilinearinterpolation is

$\begin{matrix}\begin{matrix}{{v( {x_{n},y_{n},z_{n}} )} = {{v( {x_{0},y_{0},z_{0}} )} \cdot ( {1 - {{dx}_{0}/{rx}}} ) \cdot ( {1 - {{dy}_{0}/{ry}}} ) \cdot}} \\{( {1 - {d\;{z_{0}/r}\; z}} ) + {v{( {x_{1},y_{0},z_{0}} ) \cdot {{dx}_{0}/{rx}} \cdot}}} \\{{( {1 - {{dy}_{0}/{ry}}} ) \cdot ( {1 - {{dz}_{0}/{rz}}} )} +} \\{{v( {x_{0},y_{1},z_{0}} )} \cdot ( {1 - {{{dx}_{0}/{rx}} \cdot {{dy}_{0}/{ry}} \cdot ( {1 - {{dz}_{0}/{rz}}} )} +} } \\{ {{v( {x_{0},y_{0},z_{1}} )} \cdot ( {1 - {{dx}_{0}/{rx}}} ) \cdot ( {1 - {{dy}_{0}/{ry}}} ) \cdot {{dz}_{0}/{rz}}} ) +} \\{ {{v( {x_{1},y_{0},z_{1}} )} \cdot {{dx}_{0}/{rx}} \cdot ( {1 - {{dy}_{0}/{ry}}} ) \cdot {{dz}_{0}/{rz}}} ) +} \\{{v( {x_{0},y_{1},z_{1}} )} \cdot ( {1 - {{{dx}_{0}/{rx}} \cdot ( {1 - {{dy}_{0}/{ry}}} ) \cdot}} } \\{( {1 - {{dz}_{0}/{rz}}} ) + {{v( {x_{1},y_{1},z_{0}} )} \cdot {{dx}_{0}/{rx}} \cdot {{dy}_{0}/{ry}} \cdot}} \\{( {1 - {{dz}_{0}/{rz}}} ) + {{v( {x_{1},y_{1},z_{1}} )} \cdot {{dx}_{0}/{rx}} \cdot {{dy}_{0}/{ry}} \cdot}} \\ {{dz}_{0}/{rz}} )\end{matrix} & (3.13)\end{matrix}$An alternative, more concise expression for the value of v(x_(n), y_(n),z_(n),) would be

$\begin{matrix}{{v( {x_{n},y_{n},z_{n}} )} = {\sum\limits_{i = 0}^{1}\;{\sum\limits_{j = 0}^{1}\;{\sum\limits_{k = 0}^{1}\;{( {1 - {{dx}_{i}/{rx}}} )( {1 - {{dy}_{j}/{ry}}} )( {1 - {{dz}_{k}/{rz}}} ){v( {x_{i},y_{j},z_{k}} )}}}}}} & (3.14)\end{matrix}$

Tricubic interpolation—Supersampling was also conducted using tricubicinterpolation which generates a new voxel value at location (x_(n),y_(n), z_(k),) based on the 64 nearest neighbors in the original space,v(x_(i), y_(i), z_(k)), where i, j, and k can have the values −1, 0, 1,2. Briefly, the new value can be expressed as

$\begin{matrix}{{v( {x_{n},y_{n},z_{n}} )} = {\sum\limits_{i = {- 1}}^{2}\;{\sum\limits_{j = {- 1}}^{2}\;{\sum\limits_{k = {- 1}}^{2}\;{{P( \frac{x_{n} - x_{i}}{x_{1} - x_{0}} )}{P( \frac{y_{n} - y_{j}}{y_{1} - y_{0}} )}{P( \frac{z_{n} - z_{k}}{z_{1} - z_{0}} )}{v( {x_{i},y_{j},z_{k}} )}}}}}} & (3.15)\end{matrix}$where P(t) is a cubic polynomial interpolation function with thefollowing specification:

$\begin{matrix}{{P(t)} = \{ \begin{matrix}{{1 - {2( t^{2} )} + {t^{3}}}:} & {0 \leq {t} \leq 1.0} \\{{4 - {8 \cdot {t}} + {5 \cdot t^{2}} - {t^{3}}}:} & {1 < {t} \leq 2.0} \\{0:} & {{t} > 2.0}\end{matrix} } & (3.16)\end{matrix}$Although, the inventors have used the methods set forth hereinabove, thepresent invention includes any method of supersampling which providessufficiently improved image resolution.

Section 4: Segmentation—Image segmentation is the process by which animage is divided into regions corresponding to objects and background.In the analysis of pulmonary nodules in radiographic images,segmentation involves the differentiation of nodule from non-noduleregions. This chapter details algorithmic solutions to the nodulesegmentation problem, both for well-circumscribed nodules and for thosewith attachments to vasculature or the pleural surface.

Basic Segmentation—The automated (and manual) analysis of pulmonarynodules is intimately tied to the image segmentation problem. In orderto take measurements of nodule size, shape, and density characteristics,one must first determine which image regions belong to the nodule andwhich do not. A first order solution to this problem can be describedusing a two-level model.

Two-Level Model—Considering the region of interest (ROI) as composed oftwo classes of voxels, those that are part of the nodule (foreground)and those that are not (background). FIG. 21 shows a single 2D image ofa nodule. The bright, high-attenuation region is the nodule, surroundedby a darker, low-attenuation region corresponding to lung parenchymaltissue.

Examining the frequency distribution of gray-level intensities in thisimage, it is found that they exhibit a bimodal distribution. A histogramof the image gray-levels (in HU) is shown in FIG. 22. The large modecentered at around −850 HU represents the low-intensity pixels in thebackground. The smaller mode centered near 0 HU corresponds to thehigher-intensity nodule region.

Consider the two-dimensional image function, ƒ(x, y). A simple two-levelmodel can be defined where foreground and background pixel classes maybe separated by intensity. A transformation r maps pixels in two one oftwo classes, 0 or 1. A simple transformation that relies on the twoclasses having complete separation by intensity can be described asfollows:

$\begin{matrix}{{r( {f( {x,y} )} )} = \{ \begin{matrix}{0:{{f( {x,y} )} < T}} \\{1:{{f( {x,y} )} \geq T}}\end{matrix} } & (4.1)\end{matrix}$

This gray-level intensity transformation is commonly known asthresholding. The challenge in separating the classes in this two-levelmodel, then, is the choice of an appropriate threshold value, T Severalmethods of threshold selection are described in the next section.

Thresholding—The goal of threshold selection is to determine a thresholdwhich optimally separates background from foreground regions. While thismay seem straightforward for a simple image (such as that shown in FIG.21), the optimality criteria may change in a variety of circumstances.Here are some examples: (a) maximize separation of two modes in abimodal histogram derived from a single image; (b) maximize separationof two modes in a bimodal histogram of an entire 3D ROI; (c) maintain aconsistent segmentation method between scans of the same nodule overtime; and (d) maintain a consistent segmentation method across nodulesin multiple subjects.

The threshold selection problem for a single 2D nodule image may besolved using histogram analysis and selection of a threshold value whichbest separates the two modes of the (expected) bimodal distribution.Many mathematical relationships (e.g., algorithms) have been proposedfor the global thresholding problem (a single threshold for the entireimage), as well as for local thresholding (thresholds depend on localimage properties in addition to pixel intensity) [24, 36]. A simpleglobal threshold selection algorithm will now be described, based on themethod of Ridler and Calvard [68].

Algorithm 4.1 (Iterative Threshold Selection)

Select a candidate threshold value T for histogram H(p)Δμ₁=Δμ₂=∞μ_(1p)=μ_(2p)=∞while ((Δμ₁>ε)∥(Δμ₂>ε))R ₁ =H(p): 0<p<TR ₂ =H(p): T≦p<∞μ_(1p)=μ₁μ_(2p)=μ₂μ₁=(ΣR ₁)/|R ₁|μ₂=(ΣR ₂)/|R ₂|T=(μ₁+μ₂)/2Δμ₁=μ₁−μ_(1p)Δμ₂=μ₂−μ_(2p)end

Algorithm 4.1 iteratively selects a threshold based on the meanintensities (μ₁, μ₂) of the two histogram regions (R₁, R₂) divided by T.At each iteration the new value of T is computed as the mean of the meanvalues of each region.T=(μ₁+μ₂)/2  (4.2)

The termination condition is that the means of each region change by nomore than ε, which is somewhat application dependent. The result ofAlgorithm 4.1 on the 2D image shown in FIG. 21 is a threshold of −476HU. The resulting segmented image is shown in FIG. 23.

This method works quite well for bimodal distributions where the modesencompass similar area and have minimal overlap. One factor affectingthe relative area in each of the modes, of course, is ROI selection. Ifthe nodule is surrounded by a large amount of lung parenchyma (i.e. theROI is large compared to the nodule itself), the high-attenuation modemay not have sufficient area to be distinguished above the contributionsof image noise and other high-attenuation, non-nodule structures. Thisbecomes particularly important when analyzing three-dimensional regions.

One of the issues in threshold selection is that lower threshold valuesproduce segmented nodules which appear larger, as the interface (margin)between nodule and lung parenchyma varies widely in attenuation. FIG. 24illustrates the effect of changing the threshold value T in 100 HU stepsfrom −800 HU (top left) to 100 HU (bottom right) when segmenting thenodule shown in FIG. 21.

In order to make accurate estimates of nodule growth, the segmentationshould be consistent across multiple scans taken at different times. Ingrowth estimation, the absolute error in size measurement is lesssignificant than the relative error, considering multiple studies of thesame nodule. This is due to the fact that growth estimates are computedbased on size ratios, not on absolute size. For this reason, a singlethreshold value should be used in segmentation of each of the scans.Important caveats, however, are that the scanner must be well-calibratedand the scan protocol (dose and resolution parameters) fixed for theobserved attenuation values in each scan to be comparable. Growthestimation is covered in detail in Section 5.

It may also be beneficial to attempt a standard thresholding policyacross nodules in different subjects. This is a more difficult problem,however, as mean nodule attenuation varies somewhat with nodule type.Non-solid nodules, known as ground glass opacities (GGOs), and lesionswith mixed solid and GGO components frequently require a much lowerthreshold for accurate segmentation.

A more general problem with global thresholding methods is that regionsof similar intensity may be spatially disparate, and not representcomponents of the same object. For example, nearby vessels or otherstructures of nodule-like attenuation may be included in athresholding-based segmentation. Image noise may also be of sufficientintensity to contribute to the output of a thresholding operation. Thesestructures should not be included in the nodule region to be consideredfor further size and shape analysis.

Connected Component Analysis—One solution to the problem of disparatehigh-intensity regions contributing to the segmented nodule volume isthe use of connected component labeling and analysis. In connectedcomponent labeling, each pixel/voxel is assigned an integer valueindicating the unique 2D/3D connected component of which it is a member[24]. These connected components include all voxels through which aconnected path can be followed. Such a path connects adjacent non-zerovoxels that share at least one vertex. In other words, a voxel v(x, y,z) belongs to the same connected component as any of its 26 neighbors inthe set N, where

${N( {x^{\prime},y^{\prime},z^{\prime}} )} = \begin{Bmatrix}{{x^{\prime} \ni {{x - 1} \leq x^{\prime} \leq {x + 1}}};} \\{{y^{\prime} \ni {{y - 1} \leq y^{\prime} \leq {y + 1}}};} \\{{z^{\prime} \ni {{z - 1} \leq z^{\prime} \leq {z + 1}}};} \\{( {x^{\prime},y^{\prime},z^{\prime}} ) \neq ( {x,y,z} )}\end{Bmatrix}$and the value of the neighboring voxel is non-zero.

The use of 3D connected component analysis allows noise and extraneousstructures to be removed from the nodule image data. The result ofselecting a single connected component is an object that is contiguousin 3D space.

The first step in connected component analysis is the connectedcomponent labeling itself. There two general algorithms that can be usedto accomplish connected component labeling in three dimensions, onerecursive and one iterative.

Recursive Connected Component Labeling—The recursive connected componentlabeling algorithm tales as input a three-dimensional image v andproduces a labeled image l, where each non-zero voxel value is the indexof the component of which it is a member. The algorithm has twosections. The outer section loops over all voxels v(x, y, z) in theinput image. When an object voxel is encountered (voxels with zero valueare considered as background), the component label L is incremented andthe recursive labeling function label( ) is called. This recursivefunction examines each of its 26 neighbors for connectivity. If aconnected voxel is found, it is given the same label in a recursive callto label( ).

A recursive algorithm is quite elegant in that it labels each componentcompletely and in turn. The primary disadvantage of such an algorithm,however is that each recursion places state on the system stack which isnot removed until the termination cases are reached and the recursivecalls removed from the stack.

Iterative Connected Component Labeling—To eliminate the memory demandsassociated with connective component labeling, it is possible toformulate an iterative algorithm, albeit at the cost of complexity andexecution time.

The iterative method traverses the image in order, in increasing x, y,and z. If connectivity is established between the current voxel v(x, y,z) and any of the 13 previously searched (v(x−1 . . . x+1, y−1 . . .y+1, z−1)(9), v(x−1 . . . x+1, y−1, z)(3), v(x−1, y, z)(1)), allpossible connections are recorded in a table. At the end of theiterative search through the image space, this table must be reconciledso that connected components originally assigned different labels thatare subsequently found to be connected are given a consistent label.

Size and Spatial Analysis—Once connected component labeling has beenperformed, more complete connected component analysis can take place. Insegmentation of pulmonary nodules, several selection criteria arecommonly used in connected component analysis. These criteria are usedto (i) select the component of greatest volume, (ii) select allcomponents greater than or equal to a specified volume threshold, and(iii) discard components within a specified distance of the ROIboundaries.

These three criteria are used in the following way. 1) The object ofgreatest volume in the ROI is typically the nodule. 2) In some cases,more than one object of high relative volume in the ROI may need to beselected. 3) Lastly, segmentation of nodules near other large structures(e.g. pleural surface) may not be the object of greatest volume. Inthese cases, the extraneous object is typically close to (or adjoining)the ROI boundary.

Thus, image thresholding and connected component analysis can be usedfor segmentation of pulmonary nodules when they are of relatively highintensity compared with the surrounding lung parenchyma. Such methodsare insufficient when nodules are attached to other structures ofsimilar intensity (e.g. blood vessels).

More advanced segmentation techniques are required for isolating thenodule volume in these cases. Such techniques have been developed inthis invention and are based on morphological filtering.

Morphological Filtering—Mathematical morphology is the study ofnon-linear image filters based on geometric constraints. This field wasfirst developed by G. Matheron and J. Serra in the late 1960's andpopularized in computer vision and image processing in the subsequentdecades [74, 75]. Morphological operations are useful in edge-detection,shape-detection, skeletonization, and other image segmentation tasks.The fundamental morphological operations, erosion and dilation are basedon comparison between small regions of an image and a structuringelement or kernel, translated across the image in much the same way asthe convolution kernel used in many other filtering operations.Morphological operations were first defined to operate on binary images,and whereas they have been extended to gray-level morphology, thecurrent discussion concerns the binary morphological transformations[74, 35].

Erosion is the set of all voxels in an image A that can include a smallkernel, B. It can be expressed asA⊖B={v(x,y,z)

B _((x,y,z)) ⊂A}  (4.3)where B(x, y, z) indicates the kernel B centered at voxel (x, y, z).

Dilation is the set of points where the translated kernel B(x, y, z)would touch the image A. This can also be expressed asA⊕B={v(x,y,z)

B _((x,y,z)) }∩A≠Ø  (4.4)

An additional operation of interest is reflection,A′={v(−x,−y,−z)|v(x,y,z)εA}  (4.5)which is geometric reflection of all voxels about the center of astructuring element or a complete image.

Two additional morphological operations are defined as each compositionof the previous two, as erosion and dilation do not form an inversepair. Opening, is the erosion of an image A by a kernel B followed bydilation of the result using the same kernel:A _(B)=(A⊖B)⊕B  (4.6)

Sometimes referred to as a “rolling-ball” filter, the result of amorphological opening operation can be viewed as the interior region ofan object covered by translation of the kernel into all areas where itremains completely enclosed. An opening operation has the effect ofbreaking small connections between regions and removing sharp edges.

The dual operation, closing, is a dilation followed by an erosion:A ^(B)=(A⊕B)⊖B  (4.7)

Closing can be viewed as an opening operation performed on thebackground (non-object region) of an image using the geometricreflection of the structuring kernel, B′. Graphical examples ofmorphological erosion, dilation, and opening are shown in FIG. 25. Votethe effect of the structuring element on each of the surface features inthe dilation example. An opening operation gives the resultant surfacethe characteristics of the structuring kernel; in this case, thecircular characteristics of the disk-shaped kernel. In addition, notethat to remove thin structures connected to a surface (e.g. vessels),the kernel must be of a size greater than that of the structure to beremoved. As suggested by this example, in the segmentation of pulmonarynodules, a morphological opening operation may be used to reduce thelikelihood that noise artifacts and small vessels be considered part ofa nodule.

Vascular Subtraction—Global methods for removing the pulmonaryvasculature, based on three-dimensional region growing, tree-traversalalgorithms, and other techniques have been described [90, 16, 79]. Thegoal of these global methods is to simplify the nodule detectionproblem. For the characterization of pulmonary nodules usinghigh-resolution data, however, local methods of vascular subtraction aremore appropriate. Since only a small subset of the lung volume isnormally acquired in a high-resolution diagnostic study, global methodstracing the vascular branches from the pulmonary artery and veins onwardare impossible. Even if high-resolution data of the entire lung volumewere available, this approach would likely be impractical. Furthermoreglobal or semi-global region growing schemes based on voxel intensityalone have the danger of removing small nodules when they exhibitvascular connections.

Methods—As an alternative to global region growing techniques, we use alocal filtering approach based on mathematical morphology. This methodis based on the model described in Section 2. The initial morphologicalprocessing in our segmentation algorithm consists of an openingoperation based on a spherical kernel. This kernel is passed over theinput data in a convolution-like filter for both the erosion anddilation steps, followed by connected component analysis to retain thenodule volume and discard vascular components that were initiallydisjoint or were disconnected via the opening operation.

One disadvantage of morphological opening in this application is that ithas a “smoothing” effect on the nodule surface, as surface details areremoved, such as fine spiculations. Consider the illustration in FIG.26. The illustration on the left depicts a pulmonary nodule with avascular attachment. In the center illustration, the translation of acircular kernel is shown, depicting a 2D morphological opening operation(“rolling-ball filter”). In the result shown on the right, notice thatwhile the vessel has been removed (the desired effect), some of thesurface features present in the original nodule have been deleted, asthey were smaller than the structuring kernel.

While the need to remove vessels from consideration is important, wewould prefer not to smooth away the very nodule surface characteristicswe hope to analyze. To compensate for this smoothing effect, we mayperform an iterative constrained dilation (morphological openingfollowed by geodesic dilation) process to “regrow” these features. Theentire morphological filtering process is as follows:

Algorithm 4.2 (Iterative Morphological Filtering)

Begin with an Initial Binary Image I

-   J=(I⊖S_(d))⊕S_(d) {Perform opening using a spherical kernel S_(d) of    diameter d} Perform connected component analysis, keeping the    component of interest-   while s>=2 {Number of useful dilations}-   J=J⊕Ss {Perform dilation using a spherical kernel of diameter s}-   J=J    I {Perform a voxel-by-voxel logical AND}-   S=s/2-   end

This technique restores the detailed surface features of the nodulewithout regrowing vessels more than L=d voxels from the surface. In thefirst iteration, the dilation with a sphere of diameter d extends thesurface at most d/2 voxels in the direction of the vessel. In eachsubsequent iteration, the surface may be extended by, at most, half thedistance of the previous iteration. The upper bound on the growth ofeach vessel, therefore, is a distance of d voxels, as can be seen fromthe following geometric series.

$\begin{matrix}{L = {{\sum\limits_{i = 0}^{\log_{2}d}\;\frac{d}{2^{i}}} \simeq d}} & (4.8)\end{matrix}$

In addition, the logical AND operation guarantees that all the featuresthat are grown were present in the initially thresholded image. FIG. 27illustrates the vascular subtraction using morphological opening andconnected component analysis, but without iterative filtering to regrowsurface features. Two 3D shaded surface models are shown. The left imageshows the result of a basic segmentation method based on thresholding.The right image shows the result of a 3D morphological opening using aspherical kernel of an appropriate size to remove the vessels connectedto the nodule (as well as others in the ROI), via connected componentanalysis. Note that the surface of the segmented nodule is significantlysmoothed when compared with the original thresholded data.

In comparison, FIG. 28 illustrates the result of vascular subtractionvia Algorithm 4.2, which adds the iterative morphological operationsneeded to regrow the nodule surface features. Note that the surfacefeatures present in the original thresholded data have been restored.

FIG. 29 shows another example of vascular subtraction via iterativemorphological filtering. In this example, the major vascular componentis removed, in addition to several minor vessels (note the small vesselsprotruding from the top surface of the threshold-segmented nodule(left). Again, the vascular subtraction was achieved without significanteffect on desirable nodule surface features.

Sensitivity Analysis—The two critical choices in the segmentationalgorithm for vascularized nodules are the gray-level threshold and thestructuring kernel diameter. Methods for threshold selection werediscussed hereinabove. In this section, we will examine the sensitivityof the morphological filtering methods used to remove vascularstructures from the thresholded nodule ROIs.

The diameter of the structuring kernel used in Algorithm 4.2 affectsboth the initial morphological opening of the scene, but also the sizesof kernels used to regrow the surface features. The most importantconsideration when choosing the kernel size is that the sizes of nodulesand vessels vary, as described in our vascularized nodule model (Section2). In particular, the cross-sectional diameter of attached vessels (λfrom our model), may vary considerably from case to case. Still,although we may choose a kernel of diameter d, such that it is likely tobe larger than most vessels, overestimation of the appropriate kernelsize may lead to overestimation of the nodule volume at those pointswhere considerably smaller vessels were attached. This relationship wasexpressed mathematically in Equation 4.8, where the upper bound on“regrowth” of a vessel is bounded by the diameter of the structuringkernel, d. A complete expression for this overestimate was derived forthe spherical vascularized nodule model in Equations 2.7, 2.9, and 2.10.Here, we will illustrate graphically for a three-dimensional syntheticmodel. In this example, a synthetic nodule was constructed having adiameter of 25 voxels, with an attached vessel of 5 voxels in diameter.Given our standard 0.25 mm isotropically resampled data, this wouldcorrespond to a nodule diameter of 6.25 mm and vessel diameter of 1.25mm, well within the typical range in this study.

Using this synthetic nodule, the diameter of the structuring elementused in Algorithm 4.2 was varied and the resultant segmentationsevaluated visually and numerically. Given the characteristics of theiterative vascular subtraction method, there are three regions in whichvalues of this parameter may fall: (i) d is too small—the vessel is notremoved; (ii) d is in an acceptable range—the vessel is correctlyremoved; and (iii) d is too large—the vessel is removed, but regrown toa significant degree.

Three examples from this experiment are shown in FIG. 30. On the left(d=4), the kernel is too small to remove the vessel in the initialmorphological opening. In the center image, the kernel is in theacceptable range for a proper segmentation (in this case, d=5). On theright, the kernel is somewhat too large, resulting in a significantamount of vessel “regrowth”.

An experiment was performed to assess the behavior of this segmentationtechnique using synthetic vascularized nodules (a sphere with vessels ofvarying diameters). These models were segmented using kernels ofincreasing size, and the degree of volume overestimation recorded. FIG.31 shows the relative volume of the segmented nodule model as a functionof vessel diameter, d_(v), and kernel diameter, d. Each of theseparameters are expressed with respect to nodule spherical diameter,d_(s). The data for each vessel diameter are shown beginning with thekernel size that is sufficiently large to exclude the vessel.

These data show that, although the degree of volume overestimationincreases with kernel size (over the minimum diameter required toexclude the vessel), the effect is still within 4% for kernels as largeas 60% of the diameter of the nodule and for vessels as large as 30% thenodule. Still, the overall goal is to identify that range of kernelsizes that leads to good nodule segmentation for a wide variety of invivo nodules.

To this end, additional experiments were performed to test thesensitivity of the vascular subtraction technique (Algorithm 4.2) tostructuring kernel diameter in 21 in vivo pulmonary nodules. For eachnodule, the diameter of the structuring kernel was varied from 1 to 20voxels. In each of the resultant segmentations, the degree to which thevessel or vessels (most nodules had more than one vascular attachment)were removed was evaluated.

FIG. 32 illustrates the results of these experiments. The graph showsthe relationship between structuring kernel diameter and the percentageof nodules correctly segmented. Recall that, for a given nodule and setof attached vessels, there are usually a range of values of d that leadto a correct segmentation.

In these results we may note that more than 80% of the in vivovascularized nodules were correctly segmented by using a fixedstructuring kernel diameter of 6 voxels (1.5 mm). The graph alsoexhibits a bimodal appearance. This is due to the fact that while themajority of vascularized nodules studied were connected to vesselssmaller than 2 mm in diameter, there was an additional group of noduleswith larger vascular components, resulting in the secondary peak atd=10. Therefore, although a static structuring kernel size may be chosenthat yields reasonably good results over a wide range of nodules, abetter approach is to choose a structuring kernel of the appropriatesize for a particular case.

Pleural Surface Removal—Consider the basic model of the pulmonary noduledescribed in Sections 2. It describes two classes of pulmonary nodulesthat exhibit connectivity to the pleural surface, those with apleural-tail and those that are juxtapleural. Nodules exhibiting apleural tail, a thin structure similar in size and density to the restof the nodule, can be treated in much the same way as juxtapleuralnodules.

Juxtapleural nodules, those that share a significant amount of theirperiphery with the pleural surface, require a different segmentationapproach. An effective segmentation of a juxtapleural nodule is one thateliminates pleural and other thoracic components from the nodule volume.One approach to the detection of juxtapleural nodules has been to use 2Dmorphological opening with a spherical kernel within the chest wall, andsubtraction from the lung volume [4].

Methods—A two stage approach is used in the segmentation of juxtapleuralnodules. First, the orientation of the pleural surface in the ROIcontaining the nodule is determined using three-dimensional momentanalysis [65, 64]. The method of moments is described in Section 6. Forthe purposes of the current discussion, this method is used to determinethe angles describing the orientation of the pleural surface at thepoint where the juxtapleural nodule is attached. Once this orientationhas been determined, a structuring kernel is generated with theappropriate size and orientation such that an opening operation can beused to detect the majority of the pleural surface and chest wall, whileexcluding the nodule. This kernel is disk-shaped and large enough sothat it will fit only within the chest wall region and not the noduleportion of the image. Three-dimensional image subtraction is then usedto remove these external structures from the original image. Lastly, theremaining pleural components not detected in the opening operation areremoved during the application of Algorithm 4.2. The complete method isdescribed as Algorithm 4.3.

Algorithm 4.3 (Pleural Surface Removal)

Begin with an Initial Binary Image I

-   Using moments, determine the orientation of the pleural surface-   Generate a disk-shape kernel D, oriented parallel to the pleural    surface-   J (I⊖D)⊕D {Perform morphological opening using D}-   K=I−J {Perform image subtraction}-   Continue with iterative morphological filtering {Algorithm 4.2}

An example of the pleural surface removal technique is shown in FIG. 33.The nodule, approximately 5.5 mm in diameter, is removed from thepleural surface prior to volumetric and shape characterization. Thesegmentation process was as follows.

The method of moments was used to compute the orientation of the pleuralsurface and attached nodule. Once the orientation of the pleural surfacewas determined, a disk-shaped kernel of 34 voxels (8.5) mm was generatedand used in a morphological opening to identify those voxels in thepleural wall not containing the nodule. The pleural surface componentwas subsequently subtracted from the image, leaving the nodule and asmall amount of pleural surface not identified by the opening operation.The resultant image was then segmented using iterative morphologicalfiltering using a spherical kernel of 5 voxels (1.25 mm) in diameter, toremove the remaining elements not belonging to the nodule volume. Anadditional example of this segmentation method is shown in FIG. 34. Asin the previous example, the nodule is automatically removed from thepleural surface.

Sensitivity Analysis—It should be noted that the success of the pleuralsurface removal algorithm is somewhat dependent on the geometry of thenodule ROI. This is due to the fact that estimation of the orientationof the nodule-pleura interface is based on moment calculations involvingthe thresholded ROI. In experiments performed on in vivo juxtapleuralnodules, Algorithm 4.3 was able to perform a correct segmentation in 72%of the cases. When a user was able to manually specify the orientationof the disk, the results improved, allowing for correct segmentation in83% of the cases. In the remainder of the cases, manual segmentation wasrequired. Overall, these results are promising in that they offer aconsistent, automated segmentation technique for the majority ofjuxtapleural nodules.

Three-Dimensional Surface Generation—In addition to thethree-dimensional voxel representation of a segmented nodule, it isuseful to have a polygonal surface model for the segmented volume.Surface area and curvature calculations can be much more reliablyperformed on a smoothed surface representation. Computation of thesemetrics will be discussed in Chapter 6.

Surface Tessellation—A closed 3D segmented region can be used togenerate a polygonal surface model. One such surface model is atriangular tessellation of the surface. The surface models in this workare generated using a variant of the “marching cubes'” algorithm,described by Lorensen and Cline [48]. The modifications made to thisalgorithm will be explained, following a short description of theoriginal implementation.

In the standard marching cubes algorithm, each 2×2×2 voxel neighborhoodis considered separately. The pattern of set voxels in this neighborhoodfalls in to one of 15 basic classes. These are illustrated in FIG. 35.Each of the 256 (2⁸) permutations of set voxels in a 2×2×2 neighborhoodcan be generated using these classes via symmetry, rotation, and logicalcomplement.

For each of these canonical neighborhoods, a simple set of polygonalsurfaces (composed of one or more triangles) is used to separate imagebackground from foreground (see FIG. 35). When these surface elementsare generated for each neighborhood in a 3D image, a simple polygonaltessellation of the surface is created. The appropriate polygonscorresponding to each 2×2×2 neighborhood may be assembled quickly usinga lookup table indexed by a key corresponding to one of the 256 possibleneighborhoods, or octants, in the image. A description of the process isgiven in Algorithm 4.4.

Algorithm 4.4 (Three-Dimensional Surface Generation)

-   for all N_(x,y,z)=v(x . . . x+1, y . . . y+1, z . . . z+1)    -   compute octant o(x,yz)    -   index polygon look-up table using o(x,y,z)    -   add triangles corresponding to o(x,y,z), offset by (x,y,z)-   end

Several modifications were made to the original marching cubesalgorithm. First, the internal thresholding operation and gradientcalculations have been removed. Segmentation decisions yielding thelocations of each surface boundary are made prior to polygonization, inthe 3D segmentation method. The surface gradient calculations needed forrendering Gouraud-shaded triangles are discarded, as the tessellation isused for subsequent geometric analysis. Surface rendering is handledseparately, by other tools in the VisionX system (υrend, υ3d).

While the polygons described in the original algorithm are based ontriangles, the effective aggregate polygons may be quadrilateral, ormore complex non-planar polygons. In the new implementation, each of theconstituent triangles is considered separately, so that the basicsurface unit is the triangle, simplifying subsequent surface analysis.

An additional, more interesting, change was made to the algorithm toavoid “holes” in the generated surface. When using the original marchingcubes algorithm, it is possible for holes to appear in the final surfacetessellation. It is possible to fix this problem by splitting several ofthe fifteen canonical classes and adding polygons to prevent thesediscontinuities [21].

One important note with regard to surface curvature estimation(described in Section 6) is that the triangles generated by thetessellation algorithm should be generated in a consistent right-orleft-handed system. In other words, the order in which the threevertices of each triangle is encoded should be used to differentiate theinside from the outside of the surface. This will simplify surfacenormal calculations later in the curvature estimation method.

Surface Filtering—Although the surface model of the segmented nodule isa much better approximation to the true nodule surface, it is stillsomewhat rough, as all angles are multiples of 45°. Algorithm 4.5 is atechnique for filtering the 3D surface representation. Each vertex inthe surface is replaced with the weighted sum of neighboring verticesand itself. The parameter a specifies the proportional weight given tothe neighboring vertices. When α=0.9, the new location of a vertexdepends 90% on the locations of neighboring vertices and 10% on itsoriginal location. If α is set to 1.0, none of the original vertexlocation contributes to the result. In addition, this smoothingalgorithm may be iteratively applied n times, to control the degree ofsmoothing.

Algorithm 4.5 (Three-Dimensional Surface Filtering) for k = 1 : n forall V_(i) ∈ V {for each vertex} S_(i)(x,y,z)

0 for all V_(j) ∈ adj(V_(i)) {for each vertex adjacent to V_(i)}S_(i)(x,y,z)

S_(i) + V_(j) end V_(i)

(1-a) V_(i) + α(S_(i)/|adj(V_(i))|) end end

The overall goal of the surface filtering process is to achieve arepresentation of the nodule surface that is more likely to representits actual form. In combination with the isotropic resampling techniquedescribed in Section 3, it helps to mitigate the partial volume problem.One method of measuring the accuracy of the surface representationsdescribed is to measure the surface area of each representation of anobject with known geometry. One such example is a sphere.

A sphere of diameter 10 mm was generated synthetically and sampled to anisotropic grid with 1-mm resolution. FIG. 36 shows shaded surfacerenderings of three surface representations of this sampled data. Fromright to left, they include: the original voxel surfaces, the surfacetessellation generated using the modified marching cubes technique (seeabove), and the result of one iteration of Algorithm 4.5 applied to theprevious representation.

Analytically, the theoretical surface area for this sphere is equal to100π mm²≈314.15 mm². The surface area calculated from these threerepresentations are 480 mm², 353.16 mm², and 314.931 mm², respectively.These correspond to surface area estimation errors of 52.8%, 12.4%, and0.249%, respectively. Clearly, the filtered surface approximation givesthe best estimate of the surface area. While the sphere may be anespecially good case for the surface filtering algorithm, the surface ofmost nodules can be assumed to be locally spherical (in a small voxelneighborhood). In that model, the smoothed surface representation is areasonably good one, and in any case, much more accurate than the voxelrepresentation. Further discussion of surface area estimation is givenin Section 6. The degree to which the surface is filtered, however, hasan effect on all subsequent surface-based metrics, including theanalysis of three-dimensional surface curvature distribution. This topicwill be revisited in Section 6.

Section 5: Volume Measurement—This section and much of the disclosurehereafter concerns volumetric measurements. The actual work and thediscussion accompanying the work refers to “growth” in volumetric size.However, the present invention is not limited to “growth.” It alsoincludes the basis for measuring and handling information (e.g., data)related to reduction in volumetric size which might occur as, forexample, by regression. Indeed, the scope of the invention covers“change” and volumetric size.

The most important step in volumetric measurement of pulmonary nodulesis segmentation of the nodule from the surrounding lung parenchyma andother structures. Segmentation methods were discussed in Section 4. Oncethe nodule region has been identified, volume measurement isstraightforward. The volume in voxels, is merely the sums of all “set”voxels υ in the region of interest (ROI).

$\begin{matrix}{{V({voxels})} = {{\sum\limits_{x = 0}^{M - 1}\;{\sum\limits_{y = 0}^{N - 1}\;{\sum\limits_{z = 0}^{L - 1}{{\upsilon( {x,y,z} )}\mspace{31mu}\upsilon}}}} = \{ 0 \middle| 1 \}}} & (5.1)\end{matrix}$

To give the volume standard units, we must consider the volume of eachvoxel.

$\begin{matrix}{{V( {mm}^{3} )} = {\sum\limits_{x = 0}^{M - 1}\;{\sum\limits_{y = 0}^{N - 1}{\sum\limits_{z = 0}^{L - 1}{{\upsilon( {x,y,z} )} \cdot \upsilon_{vol}}}}}} & (5.2)\end{matrix}$

In most CT data, the voxels are anisotropic, with much higher in-plane(x−y) than axial (z) resolution. The voxel volume is therefore:ν_(vol)(mm³)=υ_(xres)·υ_(yres)·υ_(zres)=υ_(in-plane) ²·υ_(zres)  (5.3)

In our system, we resample the nodule ROI to isotropic space usingtrilinear interpolation. The voxel volume is then simply the cube of theisotropic resolution.υ_(vol)(mm³)=υ_(iso) ³  (5.4)

Isotropic resampling, described in Section 3, has several benefits. Itbrings the data into isotropic space, where the resolution is the samein each dimension. This has the benefit of simplifying subsequentprocessing. Isotropic image filters may be used, rather than having tocompensate for the different resolution in the less-resolved axialdimension. More important, the data are now at a supersampledresolution, where each original voxel value is responsible for severalin the new image space. This helps mitigate partial volume effects inthe original data.

Exponential Growth Model—Exponential growth models have long been usedin the differentiation of benign from malignant nodules [56, 89, 86,80]. The basic model was developed on the assumption of uniform celldivision, with cell dividing in two, in every generation. Therefore, thenumber of cells (N) as a function of the number of generations (n) is

$\begin{matrix}{{N(n)} = {{\sum\limits_{i = 0}^{n - 1}\; 2^{i}} = {\frac{1 - 2^{n}}{1 - 2} = {2^{n} - 1}}}} & (5.5)\end{matrix}$This exponential growth can be characterized by a simple differentialequation.

$\begin{matrix}{\frac{\mathbb{d}V}{\mathbb{d}t} = {\lambda\; V}} & (5.6)\end{matrix}$The solution of this equation is straightforward.

$\begin{matrix}{\frac{\mathbb{d}V}{\mathbb{d}t} = {\lambda\; V}} & (5.7) \\{\frac{\mathbb{d}V}{V} = {\lambda\;{\mathbb{d}t}}} & (5.8)\end{matrix}$1n V=λt+C  (5.9)

Solving for the constant C at t=0:In V ₀=λ·0+C  (5.10)C=In V₀  (5.11)

Combining Equations 5.9 and 5.11 yieldsIn V=λt+In V ₀  (5.12)V=e ^((λt+In Vo)) =e ^((In V0)) ·e ^(λt)  (5.13)This yields the most common form of the generalized exponential modelfor tumor volume based on an estimate of initial volume:V=V₀e^(λt)  (5.14)where the exponential coefficient λ may be defined with respect to thenodule doubling time.

$\begin{matrix}{\lambda = \frac{l_{n}^{2}}{D\; T}} & (5.15)\end{matrix}$A simple derivation follows, considering V to be double V₀ (this occursat one doubling time).

$\begin{matrix}{\frac{V}{V_{0}} = {2 = {\mathbb{e}}^{\lambda\; D\; T}}} & (5.16)\end{matrix}$1n2=λDT  (5.17)

$\begin{matrix}{\lambda = \frac{\ln\; 2}{D\; T}} & (5.18)\end{matrix}$More generally, the exponential model can be expressed in terms of twovolume measurements Δt days apart:

$\begin{matrix}{\frac{V_{2}}{V_{1}} = {\mathbb{e}}^{\lambda\;\Delta\; t}} & (5.19)\end{matrix}$The doubling time, DT, of a nodule, then, can be computed using a ratioof two volume measurements and the time in days, Δt, between them. Thiscan be derived as follows:V₂=V₁e^(λΔt)  (5.20)

$\begin{matrix}{\frac{V_{2}}{V_{1}} = {\mathbb{e}}^{\lambda\;\Delta\; t}} & (5.21)\end{matrix}$1n(V ₂ /V ₁)=λΔt  (5.22)

Substituting our expression for λ (Equation 5.15), results in

$\begin{matrix}{{\ln( {V_{2}/V_{1}} )} = \frac{\ln\;{2 \cdot \Delta}\; t}{D\; T}} & (5.23) \\{{D\; T} = \frac{\ln\;{2 \cdot \Delta}\; t}{\ln( {V_{2}/V_{1}} )}} & (5.24)\end{matrix}$

Using Equation 5.24, doubling times can be computed for clinically-seennodules, given accurate measures of their volume at two times.Historically, however, the doubling time has been estimated usingmeasures of diameter on chest radiographs, or even in CT slices. In thiscase, there is an additional factor of three in the denominator(diameter varies as the cube root of volume).

$\begin{matrix}{{D\; T_{2D}} = \frac{\ln\;{2 \cdot \Delta}\; t}{3\;{\ln( {D_{2}/D_{1}} )}}} & (5.25)\end{matrix}$

To put the model in perspective, FIG. 37 shows a comparison ofcharacteristically malignant and benign doubling times as measured bythe increase in nodule diameter from a baseline value of 1 mm. The timescale of the graph is 10 years. It can be seen that the progression ofeach nodule from the (currently) pre-detectable size of 1 mm, variessignificantly for the three doubling times, 75, 150, and 500 days. Atthe end of the 10 year period, the benign lesion is still onlyapproximately 5 mm in diameter. The malignant lesion, however, reachesthe T1-T2 size boundary (30 mm) after only 6 years. Furthermore, it canbe seen that the malignant nodule remains in the small, sub-centimetercategory for only 4 years (≈1500 days). It is worth noting that moreaggressive malignancies, such as those with a doubling time of 75 days,will remain in the sub-centimeter category for only 2 years, and willreach 30 mm in diameter only one year later.

Gompertzian Growth Model—Although the standard model for tumor growth isexponential, other models have been proposed. These models attempt tocapture a possible retardation of growth rate for tumors of an advancedsize, and have included logarithmic [25], logistic [57], and Gompertzian[45] mathematics as a means to model this effect.

The latter model, based on a Gompertzian function, in which exponentialgrowth is retarded at larger tumor sizes, is the most common of thealternatives to the standard exponential growth model. This slowing ofthe growth rate has been hypothesized to correspond to cell death andcentral tumor necrosis, lack of nutritional supply, cell migration anddifferentiation, as well as homologous inhibition [1].

It is arguable how early in the history of lung cancer the exponentialretardation is measurable, if present at all. It has been debated in theliterature whether a Gompertzian model is clinically practical [1], orobservable in cancers outside those in animal models [89].

Much of the literature on lung cancer screening and growth analysis hasbeen based on a simple exponential growth model. In addition, the focusof our study is small pulmonary nodules, which are likely at a stageprior to appreciable growth retardation; it has been reported thatreductions from a constant growth rate may not be observable until thelesion is 2-3 cm in diameter [55]. These facts have lead to the choiceof a simple exponential growth model for our characterization of nodulegrowth.

Interscan Intervals—One of the issues surrounding doubling timeestimation is the choice of interscan interval; the number of daysbetween CT examinations. The optimization of interscan interval, t_(i)is concerned with the following two constraints: (a) t_(i) must beenough between scans so that accurate size change (growth) can bemeasured; and (b) t_(i) must not be so long to allow substantial growthof a malignancy.

Based on our experimental results regarding the accuracy of volumemeasurement as a function of nodule diameter d, (see Section 5), we candevelop an expression for optimal values of the interscan interval. Thisleads to a scheduling function for early-repeat CT (ERCT) in the studyof small pulmonary nodules.

Our data on the reproducibility of volume measurements on syntheticnodules lead to the following assumptions. For small nodules (3≦d<6),the percent error in volume estimation was 1.1% (RMS) and 2.8%(maximum). For larger nodules (6<d<11), the RMS and maximum percenterrors were 0.5% and 0.9% respectively. If we take the maximum percenterror, ε for each size and conservatively estimate that the minimumreliably-detectable percent volume change, α is twice that (2·ε), and,further, assume a conservative linear relationship between magnitude oferror and nodule size (intuitively it should be cubic, decaying muchfaster with nodule diameter), we may construct the followingrelationship α as a function of nodule diameter, d. Let ε_(S) and ε_(L)be the maximum percent errors in volume estimation for a small (d_(s))and large (d_(L)) diameter value, respectively. We assume a linearmodel, and that the minimum reliably-detectable percent volume change,α, for a given diameter, d is 2E_(d). This leads to the following modelfor α:α=md+b  (5.29)where the slope and intercept are

$\begin{matrix}{m = {\frac{{2ɛ_{L}} - {2\; ɛ_{S}}}{d_{L} - d_{S}} = \frac{2( {ɛ_{L} - ɛ_{S}} )}{d_{L} - d_{S}}}} & (5.30) \\{{and}{b = {{{2ɛ_{S}} - {m\; d\; s}} = \frac{2( {{ɛ_{S}d_{L}} - {ɛ_{L}d_{S}}} )}{d_{L} - d_{S}}}}} & (5.31)\end{matrix}$The complete expression for α, then, is

$\begin{matrix}{\alpha = \frac{2( {{d( {ɛ_{L} - ɛ_{S}} )} + ( {{ɛ_{S}d_{L}} - {ɛ_{L}d_{S}}} )} )}{d_{L} - d_{S}}} & (5.32)\end{matrix}$

For a nodule of a given diameter, d, and doubling time, DT, the minimumreliably-detectable percent volume change is, α. given by Equation 5.32.We can then reformulate Equations 5.19 and 5.15 to yield the number ofdays needed to observe a percent change in volume of α.

$\begin{matrix}{\frac{V_{2}}{V_{1}} = {{1 + ( {\alpha/100} )} = {\mathbb{e}}^{\lambda\;\Delta\; t}}} & (5.33)\end{matrix}$1+(α/100)=e ⁽1n2/DT)Δt  (5.34)1n(1+(α/100))=(1n2/DT)Δt  (5.35)Δt=DT·1n(1+(α/100))/1n2  (5.36)This interval (Δt), can be used as the interscan interval for nodules ofdoubling time DT. It is the minimum duration in days needed to observe areliably measurable volume change of α percent, given a doubling time,DT. All that remains is the choice of doubling time on which to base ourmodel.

Doubling times for most malignant nodules range from 30-400 days [47].If we are conservative, we may estimate an upper bound on the doublingtime of most malignancies as 500 days. We select a doubling time betweenthat of the of the slowest growing nodule for which we would like topredict a status of malignancy and that of the fastest growing nodulethat we would like to classify as benign. If we tailor our interscaninterval to this doubling time, nodules with shorter doubling times(more aggressive) will be detected more easily, as they will have grownproportionally more in the time between CT examinations. Nodules withlonger doubling times will exhibit little or no detectable growth inthat time and will, therefore, be characterized as benign with respectto growth rate.

Let DT_(D) be the doubling time value that separates benign frommalignant nodules. A conservative estimate of this value might be 500days. We now have all of the mathematics in place to formally state themodel for determination of the interscan interval.

Given a nodule seen at first HRCT to have a diameter of d, we maycalculate the value for the minimum reliably-detectable percent volumechange, α, using Equation 5.32, and our estimates of ε_(S) and ε_(L) atdiameters d_(s) and d_(L) for the system under study. Following thiscalculation, the interscan interval (in days) for our chosen value ofDT_(D) is given by:t _(i) =DT _(D)·1n(1+(α/100))/1n2  (5.37)Consider the following examples, based on the error estimates determinedhereinabove, and a DT_(D) value of 500 days. From our synthetic nodulestudies, ε_(S)≈3.0%, ε_(L)≈1.0%, at d_(s)=3.0 mm and d_(L)=11.0 mm.

EXAMPLE 1

A 3.5-mm (diameter) nodule is detected using HRCT. Using Equation 5.32,we determine a to be 5.75%. Equation 5.37 then yields a value of 40 daysfor t_(i).

EXAMPLE 2

A 9.0-mm (diameter) nodule is detected using HRCT. The correspondingvalue for α is 3.0%. The appropriate interscan interval is 21 days.

Each of these interscan intervals less than 90 days, perhaps the maximumtime a second high-resolution CT study should be delayed. Intuitively,one must wait a longer interscan interval to accurately assess growth ina smaller nodule. As mentioned at the beginning of this section, one ofthe constraints of this model should be to avoid the possibility ofsignificant nodule growth during the interscan interval. The maximumamount of time a small pulmonary nodule (3≧d<10 mm) will be held betweenscans, given the error estimates described is approximately 42 days,which is less than 1.5 doubling times for even the most aggressivemalignancy. Traditional exam intervals for the study of lung cancerranged from six months to two years (approximately 4-17 times as long),even though the smallest detectable lesion was at least 10 mm indiameter.

This model for interscan interval calculation can be used as thescheduling function for duration to ERCT following an initialhigh-resolution study of a pulmonary nodule detected in a lung cancerscreening program. While it may be administratively impossible toschedule patient visits on exactly the appropriate day suggested by themodel, the value of t_(i) serves as a lower bound on the time needed towait for an accurate growth assessment to be made.

Accuracy of Doubling Time Estimates—An estimate the sensitivity ofvolumetric doubling time estimation can be made using ε, the RMS percenterror in volume measurement as determined experimentally (see FIG. 38).For this model, as compared with the parameter α in the interscaninterval model, it is sufficient to use the RMS percent error ratherthan the maximum error, as we are estimating the underlying variance ofthe measurements, rather than a conservative maximum error bound. Thevolumetric doubling time calculation is a non-linear function of twovolumetric measurements which we will assume to be uncorrelated. Inaddition, it depends on the difference in time between scans, Δt, whichmay introduce additional error.

$\begin{matrix}{{D\; T} = {{f( {V_{1},V_{2}} )} = \frac{\ln\;{2 \cdot \Delta}\; t}{\ln( {V_{2}/V_{1}} )}}} & (5.38)\end{matrix}$The error in a nonlinear function, y=ƒ(x₁, x₂, . . . , x_(n)), ofuncorrelated variables can be estimated based on the variances of eachvariable. The exact differential of y,

$\begin{matrix}{{d\; y} = {{\frac{\partial f}{\partial x_{1}}d\; x_{1}} + {\frac{\partial f}{\partial x_{2}}d\; x_{2}} + \ldots + {\frac{\partial f}{\partial x_{n}}d\; x_{n}}}} & (5.39)\end{matrix}$leads to the argument that the variance in the function may be estimatedby a sum of the contributions of the variances of each variable. This isknown as the law of error propagation. The variance in the nonlinearfunction above (if all x_(i) are uncorrelated) is

$\begin{matrix}{\sigma_{y}^{2} = {{( \frac{\partial f}{\partial x_{1}} )^{2}\sigma_{x_{1}}^{2}} + {( \frac{\partial f}{\partial x_{2}} )^{2}\sigma_{x_{2}}^{2}} + \ldots + {( \frac{\partial f}{\partial x_{n}} )^{2}\sigma_{x_{n}}^{2}}}} & (5.40)\end{matrix}$In other words, the change (variance) in the function due to measurementerror (as given by the exact differential) can be expressed as afunction of the variance in measurement of each variable. We may,therefore, use this method to determine the sensitivity of the doublingtime estimation, given estimates of the variance in each measurement,V₁, V₂, and Δt.

Starting from Equation 5.38, we form the exact differential of the DTfunction. The partial derivatives of DT with respect to each variableare

$\begin{matrix}{\frac{\partial( {D\; T} )}{\partial( {\Delta\; t} )} = \frac{\ln\; 2}{\;{\ln( {V_{2}/V_{1}} )}}} & (5.41) \\{\frac{\partial( {D\; T} )}{\partial( V_{1} )} = \frac{\ln\;{2 \cdot \Delta}\; t}{{V_{1}( {\ln( {V_{2}/V_{1}} )} )}^{2}}} & (5.42) \\{\frac{\partial( {D\; T} )}{\partial( V_{2} )} = \frac{{- \ln}\;{2 \cdot \Delta}\; t}{{V_{2}( {\ln( {V_{2}/V_{1}} )} )}^{2}}} & (5.43)\end{matrix}$and the exact differential of DT is therefore

$\begin{matrix}{{d( {D\; T} )} = {{\frac{\ln\; 2}{\;{\ln( {V_{2}/V_{1}} )}}{d( {\Delta\; t} )}} + {\frac{\ln\;{2 \cdot \Delta}\; t}{{V_{1}( {\ln( {V_{2}/V_{1}} )} )}^{2}}{d( V_{1} )}} + {\frac{{- \ln}\;{2 \cdot \Delta}\; t}{{V_{2}( {\ln( {V_{2}/V_{1}} )} )}^{2}}{d( V_{2} )}}}} & (5.44)\end{matrix}$Then, using the law of error propagation, the variance in doubling tuneestimation is

$\begin{matrix}{\sigma_{D\; T}^{2} = {{( \frac{\ln\; 2}{\;{\ln( {V_{2}/V_{1}} )}} )^{2}\sigma_{\Delta\; t}^{2}} + {( \frac{\ln\;{2 \cdot \Delta}\; t}{ {V_{1}( {\ln\;{V_{2}/V_{1}}} )} )^{2}} )^{2}\sigma_{V_{1}}^{2}} + {( \frac{{- \ln}\;{2 \cdot \Delta}\; t}{{V_{2}( {\ln( {V_{2}/V_{1}} )} )}^{2}} )^{2}\sigma_{V_{2}}^{2}}}} & (5.45)\end{matrix}$Interestingly, the variance of DT is not the appropriate error bound inthis application. We would prefer, instead, to estimate the RMS error ofthe measurement, which is equivalent to the standard deviation, σ, orthe square root of the variance.

For n measurements of variable x, these are equal to

$\begin{matrix}{{{RMS}{\mspace{11mu}\;}{error}} = \;{\sigma = \sqrt{\frac{1}{N}{\sum\limits_{i = 0}^{N - 1}\;( {x_{i} - \mu} )^{2}}}}} & (5.46)\end{matrix}$Thus, the RMS error in doubling time estimation is the square root ofEquation 5.45, or

$\begin{matrix}{\sigma_{DT} = \sqrt{{( \frac{\ln\; 2}{\ln( {{V2}/{V1}} } )^{2}\sigma_{\Delta\; t}^{2}} + {( \frac{\ln\;{2 \cdot \Delta}\; t}{{V_{1}( {\ln( {V_{2}/V_{1}} )} )}^{2}} )^{2}\sigma_{V_{1}}^{2}} + {( \frac{{- \ln}\;{2 \cdot \Delta}\; t}{{V_{2}( {\ln( {V_{2}/V_{1}} )} )}^{2}} )^{2}\sigma_{V_{2}}^{2}}}} & (5.47)\end{matrix}$

Doubling Time Estimation in Incidence Cases—Estimates of volumetricdoubling time of a nodule can be made using two measurements of anodule, as described by Equation 5.24. It may also be possible toestimate the doubling time of a nodule without two volume measurements,if there is a prior study in which the nodule was undetectable. Once anodule is detected, prior CT examinations can be evaluatedretrospectively to determine if a nodule had been present. If a nodulewas in the pre-detectable phase on a prior CT, it may be considered an“incidence case,” one that became detectable between regularexaminations, as in a screening program. It is possible to estimate themaximum doubling time for these incidence cases using assumptions aboutminimum detectable size.

If we define β to be the minimum-detectable nodule diameter for a givenCT protocol, it is possible to determine an upper bound on the noduledoubling time given two examinations; one in which the nodule wasundetectable, and one in which precise volume measurements were made. Wemay assume that the nodule size was smaller than β at the time of thefirst scan and compute the doubling time using the second study in whichthe nodule volume was measured and the interval between screeningstudies, Δt_(s). The actual doubling time of the nodule would have to beequal to or smaller (faster growth rate) than the derived value. Theupper bound on the doubling time in an incidence case may be computed as

$\begin{matrix}{{DT}_{i} = \frac{\ln\;{2 \cdot \Delta}\; t_{s}}{\ln( {V_{2}/V_{\beta}} )}} & (5.48)\end{matrix}$where V_(β) is a minimum-detectable volume estimate based on theminimum-detectable nodule diameter, β. Given the assumption of aspherical nodule volume, V_(β) would be

$\begin{matrix}{V_{\beta} = {\frac{4}{3}{\pi( \frac{\beta}{2} )}^{3}}} & (5.49)\end{matrix}$As an example, consider an annual screening protocol in which 10-mmreconstructions are generated and for which the minimum-detectablenodule diameter is determined to be 2.5 mm.

EXAMPLE 3

After several annual screens, a nodule is detected in a patient that wasundetectable retrospectively in the previous exam. The nodule is scannedat high-resolution (1 mm slice thickness) for volumetric measurementyielding a volume estimate of 73.57 mm³. The upper bound on doublingtime for this nodule could be estimated using V_(β)=8.18 mm³ andΔt_(s),=365.25. The incidence estimated doubling time would then be

${DT}_{i} = {\frac{\ln\;{2 \cdot 365.25}}{\ln( {73.57/8.18} )} = {115.3\mspace{14mu}{days}}}$Such a doubling time would put this nodule well within the range ofmalignancy. Furthermore, as this is an upper bound on the doubling time,it is possible that the nodule is growing even more aggressively.

Estimation of the doubling time of incidence cases may be an importantalternative to waiting an appropriate interscan interval for a secondvolumetric measurement. This is particularly true in incidence casesthat are significantly larger than the minimum-detectable nodule size,as they are those most likely to be malignant.

EXAMPLE 4

To validate the volumetric doubling time estimation method and to makecomparisons with more traditional methods based on 2D measurement, astudy was performed using high-resolution CT data acquired at two timesfor 13 patients enrolled in the ELCAP lung cancer screening study [28].The final diagnostic status of their nodules was determined either bybiopsy or demonstrated lack of growth for more than two years (2YNC). Inaddition, three nodules were included with unknown, but likely benignstatus.

The 3D resampling and segmentation procedures described in this workwere used prior to automated measurement of nodule size. Volumetricmeasurements were made of each scan of each nodule. For comparison, areameasurements were made using the CT slice of maximum cross-sectionalarea, and diameter estimates computed from the major and minor principalaxes of the nodule in the same image.

Nodule doubling times were computed for each case based on two criteria:volumetric measurements (Equation 5.24) and area measurements (Equation5.25). Table 5.3 lists the diameter, volume, and area estimates for eachcase, as well as the two doubling time determinations and the finaldiagnostic status.

Volumetrically-determined doubling times for the malignant nodulesranged from 51 to 177 days, all well within the range for malignancy.Similarly, doubling times for the non-malignant nodules ranged from 318to 33700 days for growing nodules and from −242 to −3610 for thosenodules that were found to have reduced volume on repeat examination.With the exception of Case 14, all non-malignant nodules had doublingtimes consistent with benignity.

An example illustrating volumetric growth estimation is shown in FIGS.39 and 40. FIG. 39 shows 2D CT images of a small pulmonary nodule(baseline diameter≃5.5 mm) taken at baseline and on repeat HRCT 33 dayslater. The marked growth of this small nodule is clearly visible evenwith such a short interscan period. The nodule data were resampled to0.25 mm³ isotropic space and segmented using Algorithm 4.2, whichperforms morphological segmentation to remove the small vascularattachments. FIG. 40 shows surface-rendered 3D representations of thesegmented nodule at baseline and on repeat HRCT, respectively. Theoverall growth of this nodule is even more evident in the 3Drepresentations than from the original 2D CT slices. Note the asymmetricgrowth observable in the exaggerated protrusions on the top and frontnodule surfaces.

Volumetric analysis was performed on both scans of the nodule. It wasfound that the noble grew from 62.5 mm³ to 85.3 mm³ during the interscanperiod. Using the exponential nodule growth model and Equation 5.24, thenodule doubling time

TABLE 5.3 In-vivo nodule doubling times based on change in volume andarea measurements Volume Area DT d (mm³) (mm²) (days) based on FinalCase Δt (mm) t₀ t₁ t₀ t₁ Volume Area Diagnosis 1 36 6.9 106.9 135.7 36.536.6 104 9700 Malignant 2 20 9.3 239.8 313.8 65.9 74.1 51 78 Malignant 369 5.4 141.3 184.8 18.1 25.7 177 90 Malignant 4 71 6.5 265.2 466.4 32.666.9 87 46 Malignant 5 33 5.5 62.5 85.3 250.1 341.2 73 49 Malignant 6745 3.9 89.0 166.4 11.4 28.3 826 378 Benign 7 35 7.4 70.0 70.9 280.1283.4 2030 135 Benign 8 35 7.2 54.6 56.3 218.5 225.3 798 532 Benign 9 844.1 36.2 36.2 13.0 14.9 33700 288 Benign′ 10 225 4.0 41.5 37.6 12.2 11.8−1570 −2840 Benign 11 61 7.1 208.6 219.3 38.9 46.3 846 164 Benign 12 708.4 207.9 222.2 52.4 53.6 731 1520 2YNC 13 306 5.8 91.5 156.2 25.6 34.1396 494 2YNC 14 128 4.2 49.6 65.6 14.0 17.5 318 265 Unknown 15 140 11.9507.8 494.3 109.8 106.1 −3610 −1890 Unknown 16 111 4.6 36.7 26.7 146.9106.9 −242 −161 Unknownwas computed to be 74 days, well within the range for malignancy. Thenodule was subsequently resected and proven malignant by histologicexamination.

The asymmetry observed in this example was relatively minor. Non-uniformgrowth in three dimensions may, however, be a significant factor indoubling time estimation. When comparing volumetrically-determineddoubling times with those based on 2D area measurements, we seediscrepancies is several cases. Consider Case 1, where thevolumetrically-determined doubling time is markedly different from thatbased on area measurement. FIG. 41 shows 2D images of the CT images withmaximum cross-sectional area used for computation of area-based doublingtime.

The segmented images of the nodule at baseline and at repeat HRCT arequite similar. Their cross-sectional areas are 36.5 and 36.6 mm²,respectively, and their perimeters are 22.7 and 23.4 mm, respectively.From this two-dimensional perspective, it is impossible to discern thatthe nodule is growing in any significant manner. FIG. 42, however,illustrates the nodule in 3D, viewed perpendicular to the scanner axis.This reveals that the nodule is, in fact, growing perpendicular to thescan plane. Thus, while the area-determined doubling time for thisnodule was a benign-appearing 9700 days, the truevolumetrically-determined doubling time was computed to be 104 days,consistent with malignancy. Again, the nodule was resected and confirmedmalignant by histology. Similar cases in which 2D measurements wereinsufficient to correctly characterize growth were Cases 7, 9, and 11.True volumetrically-determined doubling times appear, therefore, to bebetter quantifiers of nodule growth, as would be expected.

Volumetric Growth Index—When evaluating the statistical differencesbetween doubling times computed for malignant and non-malignant nodules,it becomes apparent that the nonlinearity of the measure poses achallenge to evaluation. Student's t-test was used to estimate thedifference in means between the values of DT computed for malignant andnon-malignant nodules. The fact that both large doubling times(corresponding to slow growth) and negative doubling times(corresponding to reduction in nodule size) both correspond tonon-malignant processes, assessment of the means of the DT distributionsdoes not effectively characterize their differences.

To eliminate this non-linearity and to provide another intuitive measureof growth, we may define the volumetric growth index, or VGI. The VGI isdefined to be the relative increase in nodule volume in one year. Thiseffectively remaps the non-malignant nodules that grow slowly or thatare seen to be reduced in volume to the range of small positive andnegative numbers. This transformation can be accomplished as follows.Beginning with Equation 5.24, substituting one year (365.25 days) as Δt,and solving for the ratio of volume measurements we get

$\begin{matrix}{\frac{V_{2}}{V_{1}} = {\exp( \frac{\ln\;{2 \cdot 365.25}}{DT} )}} & (5.50)\end{matrix}$Subtracting one from this ratio, to give the relative increase in volumeyields the expression for VGI:

$\begin{matrix}{{VGI} = {{\exp( \frac{\ln\;{2 \cdot 365.25}}{DT} )} - 1}} & (5.51)\end{matrix}$This metric describes the fractional relative increase in nodule volumeeach year, given in units of relative volume (percent/100). Consider thefollowing three examples.

EXAMPLE 5

A benign nodule is determined to have a volumetric doubling time, DT, of1240 days. Using Equation 5.51, we determine the volumetric growthindex, VGI, to be 0.227, or 22.7%. This indicates that at the currentgrowth rate, the nodule will increase 22.7% in volume per year.

EXAMPLE 6

A malignant nodule is determined to have DT value of 75 days. We maythen determine its volumetric growth index, VGI, to be 28.2, or 2820%.In each year, this aggressive nodule with increase its volumeapproximately 28-fold.

EXAMPLE 7

A benign nodule is determined to have DT value of −2100 days. Thecorresponding value for VGI would be −0.114, or −11.4%. The nodule istherefore slowly decreasing in size.

EXAMPLE 8

With the volumetric growth index, we may now perform a statisticalcomparison of the distribution of this metric in malignant andnon-malignant nodules. Table 5.4 lists the doubling time and volumetricgrowth index for each of the in-vivo nodules given in Table 5.3.

Statistical analysis reveals that the distributions of VGI values formalignant and non-malignant nodules are significantly different. The VGIfor malignant nodules was 40.44±25.53 (mean±sem), whereas fornon-malignant nodules it was 0.2620±0.1520. The distributions showed asignificant difference in means, using Student's t-test (p<0.028). Giventhe limited sample size, and to remove assumptions about the normalityof the distributions, we may also use non-parametric statistics tocompare values for the two classes of nodules. Using Wilcoxon rank sums,the distributions were also shown to be significantly different(p<0.002).

TABLE 5.4 In-vivo volumetric doubling times and associated volumetricgrowth indices Case DT VGI I Status 1 104.4 10.31 Malignant 2 51.09140.9 Malignant 3 177.4 3.167 Malignant 4 87.04 17.33 Malignant 5 73.4230.45 Malignant 6 825.5 0.3589 Benign 7 2026 0.1331 Benign 8 797.90.3734 Benign 9 33670 0.007548 Benign 10 −1571 −0.1488 Benign 11 845.90.3489 Benign 12 730.7 0.4141 2YNC 13 395.6 0.8941 2YNC 14 317.8 1.218Unknown 15 −3606 −0.06780 Unknown 16 −241.6 −0.6493 Unknown

Section 6: Analysis of Three-Dimensional Moments—Three-dimensionalmoments [65] have been used to characterize shape information in avariety of three-dimensional data. The standard 3D Cartesian moment setof order n, m_((pqr)), contains all moments m_((pqr)), such thatp+q+r≦n. These moments are defined as follows:

$\begin{matrix}{m_{pqr} = {\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{x^{p}y^{q}z^{r}{f( {x,y,z} )}\ {\mathbb{d}x}\ {\mathbb{d}y}\ {\mathbb{d}z}}}}}} & (6.1)\end{matrix}$where f(x, y, z) is a continuous function of three dimensions. In 3Dsampled images (such as those arising in from CT data) this results in adiscrete image v (x, y, z) of size (M×N×L). The corresponding moment setcan then be calculated using the discrete expression

$\begin{matrix}{m_{pqr} = {\sum\limits_{x = 0}^{M - 1}\;{\sum\limits_{y = 0}^{N - 1}\;{\sum\limits_{z = 0}^{L - 1}\;{x^{p}y^{q}z^{r}{v( {x,y,z} )}}}}}} & (6.2)\end{matrix}$

Many of the lower order moments have physical interpretations. Thezeroth-order moment in three dimensions, for example, is the volume ofthe object:

$\begin{matrix}{m_{000} = {\sum\limits_{x = 0}^{M - 1}\;{\sum\limits_{y = 0}^{N - 1}\;{\sum\limits_{z = 0}^{L - 1}{v( {x,y,z} )}}}}} & (6.3)\end{matrix}$

Similarly, the center of mass of the object is described by (x, y, z),where

$\begin{matrix}{\overset{\_}{x} = {{\frac{m_{100}}{m_{000}}\overset{\_}{y}} = {{\frac{m_{010}}{m_{000}}\overset{\_}{z}} = \frac{m_{001}}{m_{000}}}}} & (6.4)\end{matrix}$

Geometric and Densitometric Moments—Moment analysis methods for thecharacterization of objects in images fall into two general classes,based on type of application: applications that are concerned purelywith shape characteristics, and those that take into consideration thedistribution of intensity within regions. These lead to two types ofmoment sets that can be computed as follows: (a) geometric moments—eachvoxel value is considered as binary-valued, contributing a 0 or 1 to themoment sum; and (b) densitometric moments—each voxel value is consideredto have a range of intensities, or densities.

Geometric moments are computed from a segmented, binary representationof an object or scene, whereas densitometric moments are computed from agray-level representation. Consider Equation 6.2. When computinggeometric moments, each value of v (x, y, z) is binary-valued. Incontrast, when computing densitometric moments, v(x, y, z) may have arange of values in a suitable representation (8-bit/16-bit integer,16-bit/32-bit floating point).

Both geometric and densitometric moment analyses can be used in thecharacterization of pulmonary nodules. These techniques will bediscussed in Sections set forth below.

Moment Invariants—In an effort to develop metrics that do not changewith object scale, position, or orientation, certain nonlinearcombinations of moments have been devised known as moment invariants.These are of obvious use in both pattern recognition and machine vision,as they allow objects to be described in a unique manner. One drawbackof conventional moment invariants is that they have limited physicalmeaning (as do the higher-order moments themselves).

Invariant Moment Sets—An alternative technique in object description isto establish a standard set of moments that can uniquely describe anobject. It can be shown that such a unique set would need to be ofinfinite order, but that moment sets of order 3-5 are sufficient forpractical classification of most three-dimensional objects [65, 60].

The first step toward generating an invariant moment set is to requirethe center of mass to be located at the origin.μ₁₀₀=μ₀₁₀=μ₀₀₁=0  (6.5)

This produces a central moment set, {μ_(pqr)}. If further normalizingrestrictions are placed on a moment set {m_(pqr)}, it can be transformedinto a standard moment set, {M_(pqr)}. Briefly, these restrictions areas follows:

1. Volume—The object volume must be normalized to 1.0.M₀₀₀=1.0  (6.6)

2. Center of Mass—The center of mass of the object must be at theorigin.M₁₀₀=M₀₁₀=M₀₀₁=0  (6.7)

3. Principal Axes—The principal axes of the object's ellipsoid ofinertia (EOI) must be aligned with the Cartesian axes.M₁₁₀=M₁₀₁=M₀₁₁=0  (6.8)

4. Moments of Inertia—The moments of inertia must be (in descendingorder) z, y, and x.m₂₀₀≧m₀₂₀≧m₀₀₂  (6.9)

5. Rotational Normalization—There must be a unique orientation of theobject.m₃₀₀≧0, m₀₃₀≧0, and (optionally) m₀₃₀≧0  (6.10)

Transformations of Three-Dimensional Moments—Consider the following3-by-3 matrix, specifying a generalized orthogonal transformation of the3D image space.

$\begin{matrix}{U = \begin{pmatrix}\mu_{00} & \mu_{01} & \mu_{02} \\\mu_{10} & \mu_{11} & \mu_{12} \\\mu_{20} & \mu_{21} & \mu_{22}\end{pmatrix}} & (6.11)\end{matrix}$

The effect of this arbitrary orthogonal image transformation is totransform the moment m_(pqr) to m′_(pqr). Its effect can be achieveddirectly in moment space using the following formula:

$\begin{matrix}{m_{pqr}^{\prime} = {\sum\limits_{s_{1} = 0}^{p}\;{\sum\limits_{t_{1} = 0}^{s_{1}}\;{\sum\limits_{s_{2} = 0}^{q}\;{\sum\limits_{t_{2} = 0}^{s_{2}}\;{\sum\limits_{s_{3} = 0}^{r}\;{\sum\limits_{t_{3} = 0}^{s_{3}}\;{C\mspace{11mu} U\mspace{11mu} m_{abc}}}}}}}}} & (6.12)\end{matrix}$where the combinatorial product C is defined as

$\begin{matrix}{C = {\begin{pmatrix}p \\s_{1}\end{pmatrix}\begin{pmatrix}s_{1} \\t_{1}\end{pmatrix}\begin{pmatrix}q \\s_{2}\end{pmatrix}\begin{pmatrix}s_{2} \\t_{2}\end{pmatrix}\begin{pmatrix}r \\s_{3}\end{pmatrix}\begin{pmatrix}s_{3} \\t_{3}\end{pmatrix}}} & (6.13)\end{matrix}$the transformation product U is defined asU=u₁₁ ^(t1)u₂₁ ^(t2)u₁₂ ^(s1−t)1u₂₂ ^(s2−t)2u₁₃ ^(p−s)1u₂₃ ^(q−s)2u₃₁^(t3)u₃₂ ^(s3−t)3u₃₃ ^(r−s)3  (6.14)and the indices a, b, and c of the independent moment term are asfollows.a=t ₁ +t ₂ +t ₃  (6.15)b=s ₁ +s ₂ +s ₃ −t ₁ −t ₂ −t ₃  (6.16)c=p+q+r−s ₁ −s ₂ −s ₃  (6.17)

This formulation allows scaling and rotation to be applied in one step.A complete affine transformation (translation, rotation, and scaling)and may be accomplished by translation of the center of mass to theorigin followed by a linear transformation combining scaling androtation.

Object Orientation—In order to put the object under study into standardform (and thereby compute {M_(pqr)}), it is necessary to determine itsorientation in three-dimensional space. This can be accomplished via thesolution to the following eigenproblem:Ax=λx  (6.18)where

$\begin{matrix}{A = \begin{pmatrix}m_{200} & m_{110} & m_{101} \\m_{110} & m_{020} & m_{011} \\m_{101} & m_{011} & m_{002}\end{pmatrix}} & (6.19)\end{matrix}$

The orthonormal basis of eigenvectors produced in the solution of thisproblem will point in the directions of each of the principal axes ofthe object. Using these unit vectors V_(x), V_(y), V_(z), it is possibleto find the orientation of the object. Consider a standard orientationwhere the major principal axis (V_(x),) is aligned with the x-axis, theintermediate principal axis (V_(y)) is aligned with the y-axis, and theminor principal axis (V_(z)) is aligned with the z-axis (FIG. 43).

Note that the three principal axes must be mutually orthogonal, as theyform an orthonormal basis for the object distribution. This standardorientation can also be specified in terms of second order moments(moments of inertia):m₂₀₀≧m₀₂₀≧m₀₀₂  (6.20)

The moment of inertia about the x-axis (m₂₀₀) should be greater than orequal to that about the y-axis (m₀₂₀), which should, in turn, be greaterthan that about the x-axis (m₂₀₀).

The lengths of the principal axes (twice the length of the semi-axes ofthe ellipsoid of inertia) can be derived from the eigenvalues of thesystem. If the eigenvalues of the system given in Equation 6.18 aresorted such thatλ₀≧λ₁≧λ₂  (6.21)then the lengths of the principal axes are

$\begin{matrix}{{length} = {{V_{x}} = {2{\sqrt{\lambda_{0}} \cdot {\,^{3}\sqrt{\frac{3V}{4\pi\sqrt{\lambda_{0}\lambda_{1}\lambda_{2}}}}}}}}} & (6.22) \\{{width} = {{V_{y}} = {2{\sqrt{\lambda_{1}} \cdot {\,^{3}\sqrt{\frac{3V}{4\pi\sqrt{\lambda_{0}\lambda_{1}\lambda_{2}}}}}}}}} & (6.23) \\{{height} = {{V_{z}} = {2{\sqrt{\lambda_{2}} \cdot {\,^{3}\sqrt{\frac{3V}{4\pi\sqrt{\lambda_{0}\lambda_{1}\lambda_{2}}}}}}}}} & (6.24)\end{matrix}$

Given the restrictions of standard orientation, it is possible to definenotions of roll, pitch, and yaw, allowing the specification of anarbitrary orientation as three rotations from the standard form. If werestrict objects to be longest along their major principal axis andshortest along their minor axis, these definitions correspond to theirstandard geometric counterparts. Simply, roll is a rotation about themajor principal axis; pitch is a rotation about the intermediateprincipal axis; yaw is a rotation about the minor principal axis. Usingthe principal axis vectors obtained from the solution to the aboveeigenproblem, these rotational angles are simple functions of theirprojections on the YZ, XZ, and XY planes, respectively.

$\begin{matrix}{{roll} = {{\cos^{- 1}( \frac{{proj}_{YZ}V_{z}}{V_{z}} )} \cdot {{sign}( {V_{z}(y)} )}}} & (6.25) \\{{pitch} = {{\cos^{- 1}( \frac{{proj}_{XZ}V_{x}}{V_{x}} )} \cdot {{sign}( {V_{x}(z)} )}}} & (6.26) \\{{yaw} = {{\cos^{- 1}( \frac{{proj}_{XY}V_{x}}{V_{x}} )} \cdot {{sign}( {V_{x}(y)} )}}} & (6.27)\end{matrix}$

The projections, derived from the dot products of each vector with theline of minimum distance in each corresponding plane simplify to thefollowing:proj_(yz) V _(z)=(0, V _(z)(y), V _(z))z))  (6.28)proj_(xz) V _(x)=(V _(x)(x),0,V _(x))z))  (6.29)proj_(xy) V _(x)=(V _(x)(x),V _(x)(y),0)  (6.30)

The sign ( ) (signum) function is used to uniquely define the rotations,as there are otherwise two possible solutions mapped to the same angleby the inverse cosine (the range of cos⁻¹ is [0,π]).

TABLE 6.1 Primary size metrics in two and three dimensions Feature. 2D3D gross size area volume surface perimeter surface area dimensions ofEOI length, width length, width, height

Once the orientation of an object has been determined, it is possible totransform the entire moment set {m_(pqr)} into its correspondingstandard moment set, {m_(pqr)}, using the general transformation givenin Equation 6.12.

Primary Metrics—Primary metrics of nodule size and density are thosethat are computed directly from the segmented nodule data.

Primary Size Metrics—The primary measures of nodule size in two andthree dimensions are shown in Table 6.1, where “EOI” refers to theellipse of inertia (in 2D) or ellipsoid of inertia (in 3D) representingthe nodule. Calculation of two-dimensional metrics is normally performedon the single CT image of maximum cross-sectional area.

Volume and area calculations are straightforward, as they are thezeroth-order geometric moment in 3D, and 2D, respectively. Using thegeneral equation for moments (Equation 6.2), we get volume (m₀₀₀) asexpressed in Equation 6.3. Similarly, the area of a two-dimensionalobject is m₀₀. It is important to note that these equations result insize measures in units of voxels or pixels. Section 5 discusses thenormalization of volume units based on voxel resolution.

Surface Area Calculation—There are two types of methods for surface areacalculation. Those that operate on segmented voxel data, and those thatoperate on a polygonal surface representation of the nodule (describedin Section 4 Surface area calculations made directly from the segmentedvoxel representation are basically the sum of the exposed voxel faces,those that separate object from background. The area of each exposedvoxel face is simply the product of two of the image resolution values.

Algorithm 6.1 (Surface Area Calculation (Voxel-Based)) S

0 for all v(x,y,z)

(v(x,y,z) = = 1) {for each foreground voxel} if (v(x−1,y,z) = = 0∥v(x +1,y,z) = = 0) S

S + (v_(yres) ·v_(yres)) {add yz rectangle} if (v(x,y−1,z) = = 0∥v(x,y +1,z) = = 0) S

S + (v_(xres) ·v_(zres)) {add xz rectangle} if (v(x,y,z−1) = =0∥v(x,y,z + 1) = = 0) S

S + (v_(xres) ·v_(yres)) {add xy rectangle} end

This voxel-based algorithm tends to exaggerate the surface areaestimate, as the rectangular approximation to the surface is, ingeneral, not nearly as smooth as the object it represents. FIGS. 44through 47 illustrate this point in two dimensions.

FIG. 44 shows the silhouette of an object whose shape will bediscretized by the sampling procedure. FIG. 45 shows the discretizedrepresentation of the object. Note the error introduced by sampling thecontinuous image to the grid. Similar to the partial volume problem inCT, pixels containing insufficient overlap with the object are assignedto the background. The boundary of this pixel-based representationclearly overestimates the true extent of the object surface (perimeter).FIG. 46 illustrates a polyline representation of the object boundary.The polyline (a two-dimensional analog to the polygonal surfaces used inour nodule analysis) provides a better representation of the true objectform and corresponding surface length. Finally, FIG. 47 illustrates afiltered (smoothed) polyline representation, where the location of eachvertex is replaced by a weighted sum of its original location and thatof the two adjacent vertices.

FIG. 48 illustrates the difference in voxel and surface representationsin three dimensions. Each is a shaded surface rendering of a polygonalsurface representation of the same small (d≈5 mm) pulmonary nodule. Theimage on the left is the surface given by the voxels themselves. Thecenter image is the result of surface tessellation by the modifiedmarching cubes algorithm discussed in Section 4. The image on the rightis the filtered version of the center image following two iterations ofAlgorithm 4.6.

It is clear from the figure that the voxel-based representationexaggerates the surface area. The degree to which the surface area isoverestimated can be evaluated by comparing voxel-based surface areaestimates to those computed from polygonally tessellated surfaces.

The standard surface area estimation method for this work operates onthe polygonal surface representations of the nodule described in Section4. Surface area is computed as the sum of areas of each polygon in thesurface. Algorithm 6.2 is a description of this method for surfacestessellated with triangles (the standard). Extension of the method toinclude convex planar polygons of n sides is straight-forward, as theycan be decomposed into triangles (as is the case for the squares used totessellated the voxel-based image). It is based on computing the area ofeach triangle using Heron's formula. For clarity, we define theEuclidean distance function between vertices V_(a) and V_(b) as:D(V _(a) ,V _(b))=√{square root over ((V _(a)(x)−V _(b))x))²+(V_(a)(y)−V _(b)(y))²+(V _(a)(z)−V _(b)(z)))}{square root over ((V_(a)(x)−V _(b))x))²+(V _(a)(y)−V _(b)(y))²+(V _(a)(z)−V_(b)(z)))}{square root over ((V _(a)(x)−V _(b))x))²+(V _(a)(y)−V_(b)(y))²+(V _(a)(z)−V _(b)(z)))}{square root over ((V _(a)(x)−V_(b))x))²+(V _(a)(y)−V _(b)(y))²+(V _(a)(z)−V _(b)(z)))}{square rootover ((V _(a)(x)−V _(b))x))²+(V _(a)(y)−V _(b)(y))²+(V _(a)(z)−V_(b)(z)))}²  (6.31)

Algorithm 6.2 (Surface Area Calculation (Polygon-Based)) S

0 for all T₁ = {V₀, V₁, V₂} {for each triangle} a = D(V₀, V₁) b = D(V₁,V₂) c = D(V₂, V₀) s = 1/2(a + b + c); A = {square root over (s(s − a)(s− b)(s − c))}{square root over (s(s − a)(s − b)(s − c))}{square rootover (s(s − a)(s − b)(s − c))} S

S + A end

Consider the nodule surfaces shown in FIG. 48. Using Algorithm 6.2, wefind the surface areas are 122.25 mm², 81.2759 mm², and 79.4623 m², forthe voxel-based, polygonal surface, and filtered polygonal surfacerepresentations, respectively.

This result is similar to that derived in Section 3, where it was shownthat the error in area and volume measurement is reduced with anincreasing resampling ratio. Filtering of the polygonal representationprovides much the same benefit, allowing the locations of each vertex inthe boundary to have displacements that are more finely resolved than onthe original grid. Further evidence of the benefit of surface filteringprior to the estimation of surface area was shown in Section 4. Theerror in surface area estimation of a 10 mm sphere was shown to bereduced from 12% when calculated from a polygonal surface generated bythe modified marching cubes method to 0.25% following the application ofAlgorithm 4.5.

An additional experiment was performed to assess the relationshipbetween the number of iterations of surface filtering and estimation ofsurface area. Synthetic spheres of varying diameter (from 2.5 to 25 mm)were generated, tessellated, and surface filtered a varying number oftimes prior to surface area computation. The surface area estimates werewithin 2% (mean=0.4%) of theoretical values following one filteriteration as compared with 12% (mean=9.4%) using the unfilteredtessellations. Furthermore, for spheres greater than 5 mm in diameter,the surface area estimates could be improved to within 1% (means=0.4%)following two iterations of Algorithm 4.5.

Primary Density Metrics—Computation of voxel density statistics arebased on statistical moments. More specifically, they are functions ofthe central statistical moments. Analogous to the central moments of theimage data described in Section 6.2, these moments are summations ofpowers of the voxel density values, normalized to the mean value.

$\begin{matrix}{\mu_{p} = {\sum\limits_{0}^{N - 1}\;( {{v( {x,y,z} )} - \mu} )^{p}}} & (6.32)\end{matrix}$

The number of voxels, then, is simply μ₀, the zeroth-order moment.μ₀=N  (6.33)

This is exactly analogous to area and volume in two and threedimensions, respectively. Note that the normalization of these momentsis to the mean value in the one-dimensional distribution of density, notto a particular spatial location (e.g. the center of mass). Therefore,voxel density statistics provide information about the statisticaldistribution of all voxel densities in the nodule region without regardto spatial location.

The mean voxel density is the arithmetic mean of all density values inthe segmented nodule volume and can be defined as the first-order momentdivided by zeroth-order moment.

$\begin{matrix}{\mu = {\frac{1}{N}{\sum\limits_{0}^{N - 1}\;{\upsilon( {x,y,z} )}}}} & (6.34)\end{matrix}$

Variance, a measure of the variability of density values, is derivedfrom the second-order moment, expressing the mean-squared deviation ofthe voxel values from the mean.

$\begin{matrix}{{variance} = {\frac{\mu_{2}}{\mu_{0}} = {\frac{1}{N}{\sum\limits_{0}^{N - 1}\;( {{v( {x,y,z} )} - \mu} )^{2}}}}} & (6.35)\end{matrix}$

An unbiased estimate of variance is sometimes used when it ishypothesized that the sample mean does not represent the true mean of avariable over an entire population. Such hypothesis are not usuallynecessary when dealing with the distribution of voxel values in anodule, as we often have data representing the entire nodule volume.

$\begin{matrix}{{{variance}({unbiased})} = {\frac{\mu_{2}}{\mu_{0} - 1} = {\frac{1}{N - 1}{\sum\limits_{0}^{N - 1}\;( {{v( {x,y,z} )} - \mu} )^{2}}}}} & (6.36)\end{matrix}$

The standard deviation, σ is another measure of the variation in valuesfrom the mean. It is simply defined as the square root of the variance.σ=√{square root over (variance)}  (6.37)

Two additional parameters based on statistical moments, skewness andkurtosis, are used to quantify the shape of the distribution. Skewnessmeasures the shift of the density distribution to the right or left(above or below) the mean.

$\begin{matrix}{{skewness} = {\frac{\mu_{3}}{\sigma^{3}} = \frac{\sum\limits_{0}^{N - 1}\;( {{v( {x,y,z} )} - \mu} )^{3}}{\sigma^{3}}}} & (6.38)\end{matrix}$Kurtosis is a measure of the “peskiness” of the distribution.

$\begin{matrix}{{{kurtosis} - \frac{\mu_{4}}{\sigma^{4}} - 3} = {\frac{\sum\limits_{0}^{N - 1}\;( {{v( {x,y,z} )} - \mu} )^{4}}{\sigma^{4}} - 3}} & (6.39)\end{matrix}$Note that the expression given above is normalized (by subtracting 3from the “raw” kurtosis) such that the kurtosis of a normal distributionis 0.0.

These voxel statistics can be used to assess the uniformity of densityin the nodule, which may be useful in characterizing nodule cavitation,or, more generally, cell necrosis. In addition, the mean voxel densitywithin the nodule can be used in the assessment of nodule calcification,a major predictor of benignity. If the mean density or a significantsubset of the density distribution is found to be higher than thatobserved for calcium, the likelihood of benignity is considerablyhigher. Patterns of calcification within nodules have been studied [47],but are generally only observable in nodules larger than 1 cm indiameter. In this study, however, our focus is generally restricted tosmall non-calcified nodules.

Secondary Metrics—Secondary metrics of nodule size, shape, and densityare those that are derived from primary metrics, rather than measureddirectly from the image data. One of the major benefits of secondarymetrics is that they are typically defined to be size invariant. Unlikeprimary size metrics, the size invariance of secondary metrics allowsnodules of different sizes to be compared without regard to size.Although overall nodule size is a metric typically predictive ofmalignancy, it is less appropriate as a measure of small pulmonarynodules. Size as a predictor of malignancy is typically used to stratifylarger lung masses lung masses (d>3 cm) (not generally classified asnodules) from smaller lesions. These large masses have a highprobability of malignancy, as benign nodules rarely grow to this size(82). In the study of small pulmonary nodules (d<1 cm), however, suchdifferentiation based on size is not possible. Therefore, the study ofsize-invariant metrics is concerned with determining those shape anddensity characteristics that can be used to differentiate benign frommalignant lesions when they are all of a similar size.

Aspect Ratios—The simplest secondary metrics describing the shape of anodule are the aspect ratios, simple ratios of the dimensions of thesegmented nodule volume. These dimensions are computed using momentanalysis and are (in descending order of magnitude) the length, width,and height of the ellipsoid of inertia describing the segmented volume.Computation of these dimensions is done using Equations 6.22-6.24,following the solution of the eigenproblem described above. The threeaspect ratios of these nodule dimensions are defined as

$\begin{matrix}{{LHR} = \frac{length}{height}} & (6.40) \\{{LWR} = \frac{length}{width}} & (6.41) \\{{WHR} = \frac{width}{height}} & (6.42)\end{matrix}$

Similarly, an aspect ratio may be defined in two-dimensions based on theellipse of inertia computed from a single 2D image. In this case thereis only a single ratio, that of length (major axis) to width (minoraxis).

Nodule Compactness—Another important measure of the geometricdistribution of a pulmonary nodule is compactness. Compactness is ameasure of the ratio of size to surface of a nodule.

In two dimensions, compactness has been defined as a function ofcross-sectional area and perimeter (22).

$\begin{matrix}{{Compactness}_{2D} = \frac{4{\pi \cdot {Area}}}{{Perimeter}^{2}}} & (6.43)\end{matrix}$

The constant, 47π, and exponent of perimeter are used as a normalizationsuch that the 2D compactness of a circle is equal to one.

We can extend compactness to a three-dimensional metric.Three-dimensional compactness may be defined as a function of nodulevolume and surface area [62].

$\begin{matrix}{{Compactness}_{3D} = \frac{6{\sqrt{\pi} \cdot V}}{S^{3/2}}} & (6.44)\end{matrix}$

The constant, 6√{square root over (π)}, and exponent of surface area areused to normalize the function such that the compactness of a sphere isequal to one.

Sphericity and Circularity—Two additional measures of nodule shape maybe defined that capture both the compactness and the major:minor aspectratio the nodule. Sphericity is a measure of similarity of the nodule toa sphere of the same volume.

$\begin{matrix}{{Spericity} = \frac{{Compactness}_{3D}}{LHR}} & (6.45)\end{matrix}$In two dimensions, the analogous metric is circularity.

$\begin{matrix}{{Circularity} = \frac{{Compactness}_{2D}}{LWR}} & (6.46)\end{matrix}$

Secondary Density Metrics—Secondary metrics of density distributionwithin the nodule may also be defined. The goal of these measures is toquantify the regularity of the nodule density distribution and how itdeviates from a uniformly dense sphere.

Eccentricity, ξ, of the density distribution measures the displacementbetween the geometric and densitometric centers of mass (COM). It isdefined asξ=D(COM_(geom.)COM_(dens.))  (6.47)using the D operator for Euclidean distance (Equation 6.31). To obtain atruly size-invariant metric, the eccentricity may be normalized by thecube root of the volume (an estimate of nodule radius).

$\begin{matrix}{\hat{\xi} = \frac{D( {{COM}_{{geom}.}{COM}_{{dens}.}} )}{\,^{3}\sqrt{V}}} & (6.48)\end{matrix}$

Density skew, φ_(d), is a measure of the angle between the geometric anddensitometric ellipsoids of inertia (EOI). If we define the orientationof the geometric EOI to be θ_(geom) and that of the densitometric EOI tobe θ_(dens), the density skew is simplyφ_(d)=θ_(geom.)−θ_(dens.)  (6.49)

One concern in using this metric is the stability of the measure of EOIorientation. Note that as the nodule becomes more spherical, thecalculation of the orientation becomes unstable. Therefore, the measureof density skew is less accurate as the nodule is more spherical. Thisleads to the use of the normalized density skew, φ_(d), which is simplyresealed by the sphericity of the nodule such that the less-confidentvalues are closer to zero and the more-confident values are given agreater magnitude.

$\begin{matrix}{{\hat{\phi}}_{d} = \frac{\theta_{{geom}.} - \theta_{{dens}.}}{Sphericity}} & (6.50)\end{matrix}$

Curvature Model—Analysis of surface curvature is the study of the rateof change of the surface normal, φ, with respect to the surface length.In two dimensions, this is the derivative of the normal vector withrespect to the arc length.

$\begin{matrix}{\kappa = \frac{\mathbb{d}\phi}{\mathbb{d}s}} & (6.51)\end{matrix}$

In three dimensions, an infinite number of curvatures can be defined ata point on a surface, as they are computed from the curve given by theintersection of the 3D surface with any plane containing that point andthat is parallel to the surface normal at that point.

In our discrete piecewise linear model of the 3D surface of a nodule,the situation is more straightforward. The surface curvature can becomputed as the change in surface normal between a particular vertex andany (or all) of the adjacent vertices. Curvature estimates may bedefined at each vertex, or on each triangle in the surface. Methods forthe computation of these curvature estimates will now be described.

Curvature Estimation Method—Curvature estimation methods are often basedon analysis of 2D and 3D image gradient information. In this work,curvature estimation is performed on a smoothed piecewise linear surfacemodel, derived from the results of 3D image segmentation. The surfacetessellation and smoothing are described in Section 4.

Normal and Curvature Estimation—Three-dimensional surface curvatureestimation may be performed on the nodule surface model in two steps.Curvature estimates are made for each vertex comprising the surface and,in turn, for each triangle.

For each vertex V_(i) in V, connectivity analysis is performed togenerate an adjacency list, adj(V_(i)). From this adjacency list, it ispossible to compute the set of triangles, T, of which V_(i) is a member.Consider a pair of vertices (V_(a), V_(b)) adjacent to V_(i). If theyalso share adjacency, there exists a triangle with vertices {V_(i),V_(a), V_(b)}.

This method is described more formally in Algorithm 6.3.

Algorithm 6.3 (Triangle Detection) ∀V_(i) ∈ V ∀(V_(a),V_(b))∈ adj(V_(i)){For each pair of vertices adjacent to V_(i)} if V_(b) ∈ adj(V_(a)) {Ifadj(V_(a)) contains V_(b)} T = (T∪{V_(i),V_(a),V_(b)}) end

FIG. 49 illustrates a small patch of a 3D tessellated surface. VertexV_(i), and the four vertices in adj (V_(i))={V_(a), V_(b), V_(c), V_(d)}are shown. This patch will form the basis of a running example to beused in illustrating the curvature estimation model.

Once the set of triangles, T={T_(o) . . . T_(m)}, of which V_(i) is amember has been determined, the first step is to calculate the surfacenormal for each triangle in T. If T₂={V_(i), V_(a), V_(b)}), then itssurface normal, ψ₁, can be computed as the normalized cross product oftwo of its sides:

i = V i ⁢ V a → × V i ⁢ V b →  V i ⁢ V a → × V i ⁢ V b →  ( 6.52 )

This relationship is illustrated in FIG. 50. Note that we normalize thisresult to obtain the unit surface normal vector for each triangle. Thesurface normal at vertex V₁, φ₁, can now be calculated as the average ofsurface normals of each triangle of which it is a member.

ϕ i = ( ∑ j = 0 m ⁢ ⁢ j ) /  T  ( 6.53 )This is illustrated in FIG. 51 for the vertex V. The curvature at thatvertex, V_(i), is computed using the four triangles of which it is amember, {T_(a), T_(b), T_(c), T_(d)}, and their respective surfacenormals, {ψ_(a) ψ_(b) ψ_(c) ψ_(d)}.

Curvature estimates for each vertex, C, may be derived from the vertexsurface normals, φ. The goal is to estimate the derivative of thesurface normal across points on the surface. Since the tessellatedsurface model is a piecewise linear approximation to this surface,discrete curvature estimates may be made at each vertex by consideringthe variation in vertex surface normal and that of adjacent vertices.Such an arrangement is illustrated in FIG. 52, which shows the surfacenormals at V_(i) (φ_(i)) and each of the vertices in adj (V_(i)). Notethat the surface normals of the adjacent surfaces have been determinedusing triangles not shown in the figure.

Based on this foundation of triangle and vertex surface normals, severaltypes of curvature estimates may be defined. Each of these is based onthe angular difference between the surface normal of a vertex, V_(i),and that of each of the adjacent vertices, adj(V_(i)). The angulardifference, θ_(a), between surface normals at vertices V_(i) and V_(a)can be computed using the arc cosine of the normalized scalar (dot)product of the two surface normal vectors, φ_(i) and φ_(a).

$\begin{matrix}{\phi_{a} = {\cos^{- 1}( \frac{\phi_{i} \cdot \phi_{a}}{{\phi_{i}}{\phi_{a}}} )}} & (6.54)\end{matrix}$Graphically, each of the vertices (and associated surface normals) maybe translated so that they are coincident, and the angle between thereclearly defined. This is illustrated in FIG. 53.

With this differencing mechanism in place, we may now discuss thevarying curvature estimates possible, given our 3D tessellated model.The simplest of these is the determination of the principal curvatures.The principal curvatures, κ₁ and κ₂ are defined as the maximum andminimum values, respectively, of the normal curvature evaluated at avertex. Note that in our discrete vertex space, this may be approximatedusing the extrema of the curvature values with respect to each of theadjacent vertices. Algorithm 6.4 illustrates the computation ofprincipal curvatures for each vertex in a 3D surface model.

Algorithm 6.4 (Principal Curvature Estimation) for all V_(i) ∈ V(foreach vertex) for all V_(a) ∈adj(V_(i)) θ_(a)=cos⁻¹(φ_(i)−φ_(a)) endκ₁(V_(i))=max(θ) κ₂(V_(i))=min(θ) end

Two additional curvature measures may be defined in terms of theprincipal curvatures, κ₁ and κ₂. Gaussian curvature, K, is the productof the principal curvatures at a vertex.K=κ₁κ₂  (6.55)Mean curvature, H, is the mean of the two principal curvatures.

$\begin{matrix}{H = \frac{( {\kappa_{1} + \kappa_{2}} )}{2}} & (6.56)\end{matrix}$Gaussian and mean curvatures may be computed as a simple extension toAlgorithm 6.4. While mean curvature is an excellent metric forcontinuous surfaces in 3D space, it is much more sensitive to noise inour discrete model. Therefore, a more robust measure of the averagenormal curvature at each vertex is used. Rather than computing H basedon only two atomic curvature values, we may estimate the average normalcurvature; C, as the average of angular differences in surface normalbetween a vertex and all of the adjacent vertices. This method isdescribed in Algorithm 6.5.

Algorithm 6.5 (Average Normal Curvature Estimation) for all V_(i) ε V(for each vertex) (compute curvature at vertex as average of angulardifferences) for all V_(a) ε adj(V_(j)) θ_(α) cos⁻¹((φ_(i) ·φ_(α))/(|φ_(i)||φ_(α)|)) end$C_{Vi} = {( {\sum\limits_{j = 0}^{N}\theta_{j}} )/{\theta }}$end

EXAMPLE 9

An experiment was performed to assess the improvement in accuracy of theaverage normal curvature estimation technique over mean curvatureestimation. Synthetic spheres of varying diameter were generated,tessellated using the modified marching cubes method, and filtered usingtwo iterations of Algorithm 4.5. Curvature estimates were computed basedon average normal curvature and mean curvature. The reduction in thevariance of curvature estimates was used to show the improvement incurvature estimation by average normal curvature over mean curvature, asthe true curvature of a sphere should be constant (zero variance). FIG.54 shows the results of this experiment as the percent reduction in thevariance of curvature estimates when using average normal curvature overmean curvature. There is a minimum of 50% percent reduction in variance(increase in reproducibility) when using the average normal curvatureestimates. This advantage increases with decreased sphere size, as thetessellation becomes more coarse.

Curvature estimates of each type may be made for each triangle in thesurface by averaging the curvature values for each of the constituentvertices.

$\begin{matrix}{{C_{Ti} = \frac{( {C_{Va} + C_{Vb} + C_{Vc}} )}{3}};{T_{i} = \{ {V_{a},V_{b},V_{c}} \}}} & (6.57)\end{matrix}$

This provides an elegant method for visualization of the 3D curvaturemap. Each triangle in the surface may be assigned a color (orgray-level) proportional to the curvature estimate determined for thattriangle. The complete curvature estimation technique is described inAlgorithm 6.6.

Algorithm 6.6 (Three-Dimensional Surface Curvature Estimation) for V_(i)= V_(o):V_(n) (for each vertex) for (V_(a), V_(b)) ε adj(V_(i)) if V_(b)ε adj(V_(a)) ⇒ ∃T_(j) = {V_(i), V_(a), V_(b)}(V_(i)  is  a  member  of  triangle  T_(j))end for T_(j) = T₀:T_(m) (for each triangle T_(j))$\psi_{j} = {{( {\overset{arrow}{ia} \times \overset{arrow}{ib}} )/{{\overset{arrow}{ia} \times \overset{arrow}{ib}}}}( {{calculate}\mspace{14mu}{triangle}\mspace{14mu}{surface}\mspace{14mu}{normal}} )}$end$\Phi_{i} = {{( {\sum\limits_{j = 0}^{m}\psi_{j}} )/( {m = 1} )}( {{calculate}\mspace{14mu}{vertex}\mspace{14mu}{surface}\mspace{14mu}{normal}} )}$end for V_(i) = V₀:V_(n) (for each vertex) for V_(a) ε adj(V_(i)) θ =cos⁻¹((φ_(i) · φ_(a))/(|φ_(i)||φ_(a)))|${C_{V}}_{i} = {( {\sum\limits_{j = 0}^{{N{\theta\alpha}}\;}\theta_{\alpha}} )/{N_{\theta\alpha}( {{calculate}\mspace{14mu}{curvature}\mspace{14mu}{at}\mspace{14mu}{vertex}} )}}$end end for T_(i) = T₀:T_(n) (for each triangle) T_(i) = {V_(a), V_(b),V_(c)} C_(Ti) = (C_(Va) + C_(Vb) + C_(Vb))/3 (calculate curvature ontriangle) end

FIG. 55 shows the result of 3D surface curvature estimation on amalignant small pulmonary nodule. The brighter surfaces correspond toregions of higher convexity, whereas the darker surfaces are regions ofhigher concavity. The intermediate gray-level between white and black isassigned to flat regions.

Histogram Analysis of Curvature Estimates—With methods to estimate thesurface curvature of the nodule, we may now explore how such data can beanalyzed and made suitable as a nodule metric. In particular, we wouldlike to devise a curvature metric appropriate for the differentiation ofbenign from malignant nodules.

It has been documented in the literature that malignant nodules are morebumpy (spiculated, lobular) than benign ones [26, 82, 73]. Oneappropriate goal for a curvature-based metric, therefore, would be tocharacterize the degree of surface irregularity present in pulmonarynodules. This may be accomplished using the surface curvature estimationmethod described above, combined with an analysis of the frequencydistribution of values measured in each pulmonary nodule.

Curvature Distribution Model—Consider a perfect sphere of radius r. Ateach point on its surface, the mean curvature has a constant value of1/r. The curvature of a sphere, then, is constant and inverselyproportional to its size.

In our curvature model defined in 3D tessellated space, the averagenormal curvature estimate of a “perfect sphere” asymptoticallyapproaches a single value as the density of triangles in thetessellation increases For practical spheres in our tessellated model,the average normal curvature will have a range of values whose variancedecreases with sphere size. Furthermore, the mean of this distributionwill approximate the mean curvature of the sphere, or 1/r, where r ismeasured in voxels. FIG. 56 shows the mean average normal curvatureestimate for synthetically-generated spheres of up to 10 cm in diameter,sampled on a perfectly isotropic 0.25 mm grid, tessellated usingAlgorithm 4.4, and filtered using two iterations of Algorithm 4.5. Thetheoretical curvature values by sphere diameter are shown using thedashed line in the same graph. Given the 0.25 mm resolution of the grid,the radius of each sphere in voxels is 4d/2. Consider the 1 mm sphere.Its radius is 2 voxels, and the corresponding curvature estimate is½=0.5 radians. Similarly, for the 5 mm sphere, the radius is 10 voxelsand the curvature estimate is 0.1 radians.

We may begin with a nodule model that is spherical and determine howperturbations in the nodule surface (such as those arising fromspiculations) affect the distribution of the curvature estimates made onthe nodule representation. FIG. 57 is a 2D illustration of a nodule witha surface perturbation. In this simple illustration, the surface isperturbed to contain a smaller region of constant curvature(circular/spherical). It can be seen that in this case, the frequencydistribution of curvature estimates will have three basic components, alarge peak corresponding to the general curvature of the nodule (A), andtwo additional peaks, one containing the curvature of the protrusion(B), and one corresponding to the local concavity, C, at the interfacebetween convex curvatures A and B.

Consider the following illustrative example. The surface of a sphericalnodule model of diameter d_(n) is perturbed using a local sphericalfunction at a number of points, n. The diameter of each perturbation isd_(P). Given our model, we would expect three curvature peaks in thefrequency distribution of curvature estimates for the surface, one eachcorresponding to the curvature of the original nodule surface, to thecurvature of the perturbations, and to the concavity introduced at theinterface between the large sphere and each perturbation. We willconsider an example where d_(n)=40, d_(P)=8, and n=26.

FIG. 58 shows the frequency distribution histogram) of average normalcurvature estimates in a synthetic sphere which is 40 voxels indiameter. The theoretical curvature of this sphere is therefore.1/r=2/d=0.05 radiansand the observed frequency distribution (mean=0.046, SD=0.048) agrees,as expected. Next, we may examine the distribution of measured curvatureestimates for the spheres (d_(P)=8) with which the surface will beperturbed. This histogram (mean=0.241, SD=0.032) is shown in FIG. 59.

The result of performing surface curvature analysis on the entireperturbed sphere model; however, contains more than simply the sum ofthe two characteristic curvatures of the large (nodule) and small(perturbations) spheres. FIG. 60 shows the frequency distribution ofcurvature estimates for the entire model. In addition to the peakscorresponding to the nodule surface (0.05 radians) and to the curvatureof each perturbation (0.25 radians), there is also a concave componentat approximately −0.2 radians. This is, in fact, the manifestation ofthe expected concavities at the interface between the spherical surfaceand each smaller spherical perturbation.

The location of the two peaks added to the histogram as a result ofsurface perturbation increase in absolute magnitude with decreasingperturbation diameter. The convex component increases due to itsdecreasing radius of curvature, and, correspondingly, the magnitude ofthe (negative) concave component increases as the angular differencebetween the surface and perturbation increases at the point of junction.As an example, FIG. 61 illustrates another perturbed sphere model wherethe diameter of the perturbations is half that used in the example shownin FIG. 60. Note that the curvature peaks for the perturbation surfaceand the nodule-perturbation interface have each moved away from theorigin.

Using such a model, several observations can be made. First, it can beseen that as the number of surface perturbations increases, the numberof added curvature peaks increases. Second, the curvature of each of theperturbations is higher than that of the underlying nodule surface.Third, the curvature estimates resulting from the interfaces of eachoutward perturbation contribute negative (convex) curvature values tothe distribution. Each of these observations lead to the conclusion thatwith increased numbers of surface perturbations and decreasing diameterof each perturbation, the width (and correspondingly, the variance) ofthe curvature distribution will increase. Thus, a measure of thevariance of the frequency distribution of curvature estimates willprovide a quantitative means of characterizing the irregularity of thenodule surface. Furthermore, as nodule surface irregularity has beenshown to correlate with malignancy, this curvature variance metric mayprove useful in the discrimination of benign and malignant nodules.

One more note must be made on the evaluation of the distribution ofcurvature estimates. When analyzing the 3D curvature of a nodule whichhas been segmented to remove other structures (e.g. vessels), there maybe artifacts, such as artificially convex regions, at these points ofsegmentation. For this reason, when 3D curvature analysis is performed,the curvature estimates for those regions should be discarded. Theprocedure may be done as follows: (a) segment the pulmonary nodule,removing extraneous structures; (b) generate the 3D tessellated surfacemodel; (c) compute 3D curvature estimates for each vertex and trianglein the model; (d) mask the surface model of the nodule using thedilation of each of the deleted components; and (e) compute thefrequency distribution of the remaining curvature values. This was themethod used in the 3D curvature analysis of all of the nodules in thiswork. The results are described in below.

Sensitivity of Curvature Estimates—Curvature estimates and the analysisof their distribution are sensitive to three general parameters: (i)image resolution; (ii) segmentation method; and (iii) degree of surfacesmoothing

The original CT image resolution and the isotropically-resampled imageresolution affect the degree of accuracy with which curvature estimatescan be made. Improved (finer) resolution allows curvature estimates tobe made more locally, at the cost of increased noise. Parameters of thesegmentation method, including any intensity thresholds used, affect thecurvature estimates by redefining the boundary of the nodule beingstudied. Finally, the degree of surface smoothing (in this case, thenumber of iterations of Algorithm 4.5 affects the curvature estimates byremoving high-frequency noise due to spatial quantization.Over-application of this smoothing process, however, may result in theelimination of true curvature characteristics.

EXAMPLE 10

Experiments were performed to assess the affect of Algorithm 4.5 on thedistribution of surface curvature estimates. First, a synthetic sphereof 10 mm in diameter (40 voxels) was tessellated using Algorithm 4.4. Avarying number of smoothing iterations were then applied and thevariance of the surface curvature distribution calculated. FIG. 62 showsthe results of this experiment. It can be seen that in each iteration,the variance of the curvature distribution decreases as each curvatureestimate approaches the true mean value for the sphere. This observedreduction in variance is approximately proportional to 1/n₂, where n isthe number of smoothing iterations.

Malignancy Models—Malignancy models are mathematical tools used topredict the likelihood of nodule malignancy based on a variety ofparameters. The majority of these parameters may be calculated featuremetrics, such as those described in the preceding sections. Malignancymodels may also include, however, risk-related patient demographicinformation, such as age, smoking history, and exposure to other knownlung carcinogens. The actual selection of parameters may, therefore, bebased on a number of methods, mostly statistical, that maximize theability of the model to characterize the likelihood of nodule malignancyacross a wide range of pulmonary nodules.

Statistical Parameter Selection—Parameters in a model of nodulemalignancy are generally selected in two steps. The first step that isemployed is to include parameters that would be expected to be differentin benign and malignant nodules, based on the underlying biology. Thesecond step is to select parameters through the use of statistical teststhat quantify the degree to which each parameter is associated withmalignancy.

Student's t-test was used to compare the distributions of values foreach shape parameter in malignant and non-malignant nodules. This servedtwo purposes. The first was to determine which parameters, on there own,were significantly different in the two classes of nodules. Thisinformation is quite useful in that it can lead to new research intorefinement of metrics that appear promising as predictors of malignancy.The second goal was to determine those parameters which might besuccessfully combined into a model that would be predictive ofmalignancy. Once promising individual parameters were determined,combinations thereof can be used to generate logistic regression modelsfor the classification of nodules.

The logistic regression model is a type of nonlinear multiple regressionwhich is used to model data with categorical (non-continuous outputs).For example, a standard multiple linear regression model, y=f(x₁, x₂, .. . , x_(n)) uses a linear combination of the input parameters, x_(i),to predict the output y:y=a+b ₁ x ₁ +b ₂ x ₂ + . . . b _(n) x _(n)  (6.58)

The regression method, generally based on linear least-squares methods,uses a known set, of observation vectors, x, and the correspondingvector of outputs, y, to optimize the coefficients of the model, a andb₁ . . . b_(n). The range of output values in linear regression model iscontinuous, as the output, y could be any real number.

In the case where the outputs are categorical, however, logisticregression models may be used. These allow the output to be restrictedto a particular range, typically between 0 and 1. For cases where theprocess being modeled has multiple discrete values, multinomial (orpolytomous) logistic regression may be used. However, most frequently,as in the case of nodule malignancy models, the dependent variable isbinary-valued, either malignant or non-malignant (benign). Therefore, abinary logistic regression model is appropriate.

Logistic regression is based on the logit transformation, which is thenatural log of the odds ratio. The odds ratio of an event is a functionof the probability, P of that event:

$\begin{matrix}{{odds} = \frac{P}{1 - P}} & (6.59)\end{matrix}$The logit transformation, then is

$\begin{matrix}{{{logit}(P)} = {{1{n({odds})}} = ( \frac{P}{1 - P} )}} & (6.60)\end{matrix}$

Logistic regression models are based on the assumption that thedependent variable (e.g. malignant or benign) is linearly dependent noton the observed parameters, x_(i), but on the logarithm of the oddsratio. Therefore,

$\begin{matrix}{{1{n({odds})}} = {{1{n( \frac{P}{1 - P} )}} = {a + {b_{1}x_{1}} + {b_{2}x_{2}} + {\ldots\mspace{14mu} b_{n}x_{n}}}}} & (6.61)\end{matrix}$It can be seen, then, that the probability of the event in terms of theobserved parameters is

$\begin{matrix}{P = \frac{\exp( {a + {b_{1}x_{1}} + {b_{2}x_{2}} + {\ldots\mspace{14mu} b_{n}x_{n}}} )}{1 + {\exp( {a + {b_{1}x_{1}} + {b_{2}x_{2}} + {\ldots\mspace{14mu} b_{n}x_{n}}} )}}} & (6.62)\end{matrix}$which is equivalent to

$\begin{matrix}{P = \frac{1}{1 + {\exp( {- ( {a + {b_{1}x_{1}} + {b_{2}x_{2}} + {\ldots\mspace{11mu} b_{n}x_{n}}} )} )}}} & (6.63)\end{matrix}$Thus, if our logistic model is

$\begin{matrix}{y = \frac{1}{1 + {\exp( {- ( {a + {b_{1}x_{1}} + {b_{2}x_{2}} + {\ldots\mspace{11mu} b_{n}x_{n}}} )} )}}} & (6.64)\end{matrix}$the range of outcomes is between 0 and 1 for any values of the observedparameters, and furthermore, the output of the model varies sigmoidallybetween the two extremes.

The logistic regression model was used in this work to developmalignancy models based on two- and three-dimensional shape and densityfeatures. Selection of the individual parameters and performance of themodels will be discussed in Section below.

EXAMPLE 11

In order to quantify the ability of each of the shape metrics studied incharacterizing the differences between benign and malignant smallpulmonary nodules, a series of tests were performed using a group of invivo pulmonary nodules, scanned using high-resolution CT. The metricswere first compared individually, to assess differences in thedistributions of each parameter in the two nodule classes. Next,logistic regression models were developed based on 2D and 3D metrics, aswell as metrics derived from 3D surface curvature analysis.

Individual 2D and 3D Metrics—As a first step in assessing thesuitability of each shape metric for the prediction of malignancy,Student's t-test was applied to test the difference in means of thedistribution of each parameter for benign and malignant nodules. Table6.2 shows the significance of the difference in means between values ofsize-variant 3D shape metrics for 22 benign and malignant pulmonarynodules In addition, the means and standard deviations of thedistribution of each metric for benign and malignant nodules is given.As size has been well-documented as associated with malignancy, it isnot surprising that each one of these metrics significant (p<0.05), andthat most are quite significant (p<0.01).

TABLE 6.2 Analysis of distribution of size-variant 3D metrics for benignand malignant nodules Benign Malignant Parameter P Value Mean SD Mean SDheight 0.0015 3.53413 0.81352 6.02052 2.39888 surface area 0.0031 96.30962.419 339.126 271.547 x-extent 0.0037 5.7000 1.88557 10.6429 5.26952z-extent 0.0038 5.05000 1.39578 8.64286 3.81842 width 0.0054 5.163281.63050 8.21210 3.00452 length 0.0081 6.6068 2.47830 10.7895 4.22890volume 0.0087 75.546 66.686 383.776 410.535 mass 0.0096 11931.1 12136.257878.3 61302.5 y-extent 0.0178 6.21667 2.43645 9.78571 4.06568

TABLE 6.3 Analysis of distribution of size-variant 2D metrics for benignand malignant nodules P Benign Malignant Parameter Value Mean SD Mean SDlength 0.0010 5.73953 1.94499 9.96925 3.20678 x-extent 0.0071 5.466671.83193 9.71429 4.90626 y-extent 0.0194 5.43333 2.02543 8.42857 3.53764mass 0.0240 1086.72 627.56 2357.92 1842.10 width 0.0247 4.57575 1.185936.70308 2.98820 area 0.0444 36.1917 37.4608 81.4286 61.6513 perimeter0.2400 33.7924 43.8789 66.6005 84.7562

For comparison, Table 6.3 shows the results for the same 22 nodulesusing 2D size-variant features computed on the image of maximumcross-sectional area for each nodule. Note that the distributions ofeach of the nine 3D metrics exhibit more significant differences betweenbenign and malignant nodules than all but two of the five 2D metrics.

As the focus of this work is on small pulmonary nodules, such as thosethat may be detected in a lung cancer screening program, we may notalways expect size to be a useful parameter. This is due to the factthat nodules may be detected at any point in their growth progression.Therefore, malignancies may frequently be caught while they are stillsmall (especially as incidence cases, see Section 5). Furthermore, inthose cases that are of larger size (i.e. diameter>3 cm), malignancy isless often in doubt, and biopsy is more straightforward. For thesereasons, this study is primarily restricted to nodules with diameterssmaller than 1 cm. With that in mind, size-invariant metrics shouldoffer a better hope of differentiating small benign and malignantnodules that are all of a similar size. Table 6.4 shows the results oft-test analysis of size-invariant 3D shape metrics applied to 22 smallpulmonary nodules. The metrics are listed in ascending order of p value(descending order with respect to significance).

Results for a comparative group of size-invariant 2D metrics are shownin Table 6.5. When possible, the analogous 2D metrics to those studiedin 3D were included (e.g. compactness). Here, the difference between 3Dand 2D shape metrics is less clear. As with the 3D size-invariantmetrics, there are relatively few that exhibit high significance whenused as the sole parameter in differentiating benign from malignantnodules. Therefore, we may turn to another, more appropriate measure ofthe relative benefit of 2D versus 3D nodule shape metrics. Groups ofeach may be included in a logistic regression model and the relativequalities of the models compared.

Logistic Regression Models—To test the ability of multiple featuremetrics to characterize differences between benign and malignantnodules, logistic regression models were developed based on metricswhose individual t-tests showed promising differences in distributionbetween the two nodule classes. In addition, the models were restrictedto contain parameters that captured uncorrelated information. Modelswere developed using 2D and 3D metrics to compare the ability of each ofthese groups in the characterization of nodule malignancy. Each wasbased on three features.

The 2D model included LWR, an aspect ratio based on the ellipse ofinertia of the nodule cross-section; APR, an alternate measure of 2Dnodule compactness; and voxel skewness, a measure of the normality ofthe distribution of density within the nodule.

TABLE 6.4 Analysis of distribution of size-invariant 3D metrics forbenign and malignant nodules Benign Malignant Parameter P Value Mean SDMean SID normalized volume/surface ratio 0.0053 0.42401 0.01615 0.389160.03701 3D compactness 0.0064 0.81399 0.09248 0.64124 0.17669volume/surface ratio 0.0124 0.69891 0.14654 0.96051 0.30697 sphericity0.1861 0.47675 0.16617 0.37546 0.15039 {circumflex over (φ)}_(d) 0.18886.0125 14.043 50.9941 130.101 XYR 0.1984 0.96358 0.12552 1.08233 0.29998φ_(d) 0.2234 4.0257 11.4424 25.9638 67.4143 yaw 0.2887 −70.419 41.1223−49.083 46.3801 variance 0.3676 445.517 248.491 603.033 566.038 roll0.3846 −3.1443 75.422 −38.859 111.360 ξ 0.5462 0.03285 0.02325 0.040760.03711 WHR 0.6037 1.48294 0.45815 1.38355 0.27398 pitch 0.6100 24.685118.9031 31.1248 40.2918 {circumflex over (ξ)} 0.6121 0.00906 0.008160.00716 0.00782 LWR 0.6519 1.27539 0.18937 1.31137 0.12034 YZR 0.71011.20822 0.27464 1.14632 0.50272 standard 0.7378 20.2685 6.0980 21.602812.6127 deviation mean 0.7424 144.513 40.7856 137.408 57.8951 LHR 0.74631.86727 0.52547 1.79719 0.28645 skewness 0.7777 0.07568 0.32579 0.118360.32586 kurtosis 0.9956 2.20742 0.36081 2.20840 0.44150

TABLE 6.5 Analysis of distribution of size-invariant 2D metrics forbenign and malignant nodules Benign Malignant Parameter P Value Mean SDMean SD LWR 0.0046 1.24344 0.17693 1.58450 0.32985 circularity 0.03130.64667 0.26691 0.37888 0.21562 skewness 0.0852 −0.626 0.52973 −0.15230.65853 2D compactness 0.0875 0.79975 0.30742 0.55109 0.28999 apr 0.09531.18576 0.31016 1.56414 0.72021 napr 0.1722 0.24347 0.06839 0.198260.07285 ξ 0.1994 0.04062 0.02536 0.06981 0.07872 XYR 0.2695 1.042460.222191 1.9934 0.43396 φ_(d) 0.4483 4.23139 11.3991 0.84959 0.9959{circumflex over (φ)}_(d) 0.4678 43.7186 146.571 2.1665 2.1970 yaw0.4889 −42.213 45.7468 −24.268 73.6191 mean 0.6055 152.847 43.1481141.136 59.8249 variance 0.7178 481.253 305.517 544.265 502.266 kurtosis0.7660 0.33915 25.7179 3.31157 0.8487 standard 0.9083 20.9205 6.833820.4526 12.1223 deviation {circumflex over (ξ)} 0.9466 0.00166 0.001340.00160 0.00298

The 3D model included 3D compactness, a gross measure of surfaceirregularity; normalized density skew ({circumflex over (φ)}_(d)), ameasure of the irregularity of the density distribution within thenodule volume; and the variance of the average normal surface curvaturemetric, a fine measure of surface irregularity.

The pertinent shape metrics were computed for each of 22 small pulmonarynodules (7 malignant, 15 non-malignant). Statistical analysis wasperformed to assess the ability of each model to differentiate malignantfrom non-malignant nodules. The results of this experiment aresummarized in Table 6.6.

TABLE 6.6 Performance of 2D and 3D metric-based models in characterizingmalignant and non-malignant small pulmonary nodules 2D Metrics 3DMetrics Nodule Status Malignant Benign Malignant Benign Malignant (7) 6 1 FN 7 0 Non-Malignant (15)  1 FP 14  1 FP 14 R² 0.690 0.778 P <0.0003<0.0001

The model based on 2D metrics was able to characterize much of thedifference between benign and malignant nodules. The R₂ value (used toassess the proportion of the data explained by the model) was 0.6900,and the p value (probability that the effect would be produced bychance) was 0.0003. The model was able to classify all but onemalignancy (1 false negative (FN)), and one benign nodule (1 falsepositive (FP)) correctly. Therefore, the sensitivity was 86% and thespecificity was 93%. The positive predictive value was 86%.

The model based on 3D metrics fit the experimental data even moreclosely. The value of R₂ was 0.7776, with p<0.0001. In addition, themodel was able to classify each of the cases correctly, with theexception of one benign nodule (1 FP). The sensitivity was 100%, thespecificity was 93%, and the positive predictive value was 88%. This isan especially good result considering that the penalty in producing afalse-positive result (unnecessary biopsy) is much less severe than afalse-negative one (missed cancer).

An n-way (leave-one-out) cross-validation approach was also used toevaluate the performance of the model based on 3D metrics, as it mayproduce a somewhat less-biased evaluation of the model. In thisexperiment, the model was able to correctly classify 86% of the nodulescorrectly, producing 2 FP and 1 FN. The sensitivity was 86%, thespecificity was 87%, and the positive predictive value was 75%. Toprovide a more general evaluation, Receiver-Operating Characteristic(ROC) curves were generated based on the ability of the model tocharacterize the entire dataset, as well as each of the individualmodels developed in the cross-validation scheme [27, 2]. The area underthe ROC curve, A_(z), for the overall model was 0.9905. The ROC curvefor the complete dataset is shown in FIG. 63. The mean area under theROC curves using the cross-validation scheme was 0.9870, with 95%confidence intervals between 0.9822 and 0.9918, while the area under theaggregate curve was 0.9203. This aggregate ROC curve is shown in FIG.64. Overall, the model showed considerable ability to characterize thecases in the study.

EXAMPLE 12

Analysis of 3D Surface Curvature—To evaluate the usefulness of the 3Dsurface curvature analysis itself, an experiment was performed using theset of 22 in-vivo pulmonary nodules used to evaluate the 2D and 3Dfeature metrics. The measure of the three-dimensional surface curvaturewas computed using the techniques described above, including theanalysis of surface curvature distribution. Each of the nodules studiedhad complete 3D information derived from high-resolution (1 mm slicethickness) helical CT. Curvature values computed from portions of thesurface in contact with other structures (e.g. vessels) were removedfrom consideration. Frequency distribution of this measure as well asthe mean, variance, coefficient of variation, and skewness weredetermined, as described hereinabove. The distribution of each of thesemetrics in the 22 in-vivo pulmonary nodules were examined. The resultsare summarized in Table 6.7.

TABLE 6.7 Analysis of distribution of 3D curvature-based metrics forbenign and malignant nodules Benign Malignant Parameter P Value Mean JSDMean SD variance <0.0001 0.01084 0.00395 0.02274 0.00499 standard<0.0001 0.10234 0.01983 0.15001 0.01653 deviation mean 0.0003 0.103490.02682 0.05095 0.02419 range 0.0006 0.85084 0.26098 1.94865 0.99224skewness 0.0059 0.35311 0.55021 −0.5373 0.78756 kurtosis 0.0492 4.407532.40237 6.99628 3.29447 maximum 0.0999 0.58676 0.19492 0.85648 0.54787

The mean curvature estimate in malignant nodules was 0.0509±0.010(mean±sem) vs. 0.1035±0.007 in benign nodules (p<0.0003). The varianceof curvature in malignant lesions was 0.0227±0.002 vs. 0.0108±0.001 inbenign ones (p 0.0001). The skewness of the curvature distribution was−0.537±0.298 vs. 0.353±0.142 in malignant and benign nodules,respectively (p<0.01). Finally, as a mean-normalized measure ofcurvature distribution, the coefficient of variation (CV) of curvatureestimates for malignant nodules was 3.823±0.919 vs. 1.095±0.132 inbenign nodules (p=0.0004). Thus, each metric of the measure of 3Dcurvature was useful in differentiating between small malignant andbenign nodules.

Region Analysis Without Segmentation—Image segmentation is one of themost challenging steps in the analysis of small pulmonary nodules. Itmay be difficult to remove connected structures (e.g. vessels) withoutremoving information about the surface of the nodule itself. Inaddition, when studying the change in the nodule over time (as involumetric doubling time estimation), it is of great importance that thesegmentations of each image are consistent.

It may be possible to study the change in nodule size and shape overtime, without performing explicit segmentation, however. This chapterdescribes techniques for comparing three-dimensional nodule regions ofinterest (ROIs) as two density distributions, without segmentation.

Three-Dimensional Region Registration—Prior to comparison of twothree-dimensional ROIs, they must be well-defined so that eachcorresponds to the same lung region in each study. This is a problem ofthree-dimensional region registration. Two techniques for regionregistration will be discussed: 1) center of mass (COM) alignment and 2)correlation-based region matching.

Center of Mass Alignment—The two 3D ROIs may be aligned using a simpleprocedure based on the computation of the center of mass. Beginning froman initial starting point in the 3D image, an iterative search method isused to find the location of the COM in each of a number of translatedROIs. In each iteration, the current location of the COM, C_(c) iscompared with the actual COM measured for the current ROI, C_(m). Thenext location of the COM is then replaced with C_(m) and the ROIrecomputed. The stopping condition is when the location of the COM doesnot move by more than some small distance ε between iterations.Algorithm 7.1 describes the center of mass determination technique,where D(V_(a), V_(b)) is the Euclidean distance function between voxelsV_(a) and V_(b).

Algorithm 7.1 (Iterative Center of Mass Determination) Select an initiallocation for the COM: C_(i) = v(x_(i), y_(i), z_(i)) Select the desiredsize of the ROI: S = x, y, z ΔC = ∞ C_(c)= C_(i) while (ΔC > ε) Computenew ROI bounds based on C_(c) and S C_(m) (x) = m₀₀₀/m₀₀₀ C_(m) (y) =m₀₁₀/m₀₀₀ C_(m) (z) = m₀₀₁/m₀₀₀ ΔC = D(C_(ml) C_(c)) C_(c) = C_(m) end

The inputs to this procedure, in addition to the 3D image, are aninitial starting point for the search, and the dimensions of the ROIdesired. A two-dimensional example of the operation of Algorithm 7.1 isshown in FIG. 65. Each image represents one iteration of the algorithm.In each frame, the white rectangle illustrates the current ROI, in whichthe white cross shows the computed center of mass, C_(m).

This value then becomes the center of the ROI in the next iteration. Thealgorithm stops after 10 iterations, as the COM did not change by morethan ε=2 voxels between the final two iterations.

Algorithm 7.1 can be used for the determination of two matching ROIs insequential scans of a pulmonary nodule, allowing for comparativeanalysis. The method is quite sensitive to the exact distribution ofdensity, including that corresponding to structures in the periphery(such as vessels), however. In particular, when a nodule exhibitsnon-uniform growth, the resultant ROIs may not correspond perfectly.

Correlation-Based Region Matching—In the method described in theprevious section, we aligned two images by determining a robust centerof mass for each and aligned these, using the assumption that the COMwould not shift significantly between scans. An alternative approach tothe determination and alignment of two comparable ROIs may be defined interms of image correlation. The two ROIs may be aligned using 3Dcorrelation as a measure of image correspondence. The method requiresselection of the two ROIs and specification of the limits on the extentof the search. The second ROI is then translated to all locations in thesearch region and the 3D correlation computed between the two ROIs. Thetranslation yielding the greatest value of this match metric is used tospecify the alignment of the ROIs for subsequent analysis. The techniqueis described in Algorithm 7.2, where the T (ROI, i, j, k) specifies thetranslation of an ROI coordinate system by x=i, y=j, z=k, and the star(★) operator is the three-dimensional correlation described in Equation7.1.

$\begin{matrix}{{c( {i,j,k} )} = {\sum\limits_{x}{\sum\limits_{y}{\sum\limits_{z}{{f( {x,y,z} )} \cdot {g( {{x - i},{y - j},{z - k}} )}}}}}} & (7.1)\end{matrix}$

In each iteration, the correlation of the two ROIs at the currenttranslation of ROI_(b) is computed. If the value of this correlation, c,is greater than the previous maximum value, M, the value and thelocation of this better match, L, are recorded. The final result of thealgorithm is that translation, L, that produces the best alignmentbetween the ROIs.

Algorithm 7.2 (Correlation-Based Region Matching) Select two ROIs:ROI_(a), ROI_(b) Select the extent of the search region in eachdimension: S = {x, y, z} M=0 for k = − S_(z) : S_(z) {for eachtranslation in search region} for j = −S_(y) : S_(y) for i = S_(x) :S_(x) c(i, j, k) = ROI_(a)

ROI_(b) if c(i, j, k) > M M = c(i, j, k) L = (i, j, k) end end end end

An example of the result of correlation-based registration is shown fortwo scans of a small pulmonary nodule in FIG. 60. When usingcorrelation-based matching, changes in the absolute center of mass ofthe nodule are less likely to affect the registration, as the overallmatch is based on direct correspondence between areas of intensity ineach ROI. For this reason, the correlation-based matching technique waschosen for the analysis methods described later in this chapter.

Three-Dimensional Region Weighting—Since we are analyzing the nodule ROIwithout performing explicit segmentation, vessels and other confoundingstructures will influence the metrics computed on the region as a whole.In an effort to reduce the influence of peripheral structures, we mayprefilter the region using a weighting function that gives the highestweight to central structures (near the center of mass), and less weightto the periphery. In this section we will examine several weightingfunctions and their effects on our model of the nodule region.

Weighting Functions—Two region weighting functions were considered whendeveloping the region analysis technique. These included: (a)rectangular radial window; (b)triangular radial window; and (c) Guassianwindow.

In the following discussion, we will assume the ROI to be translatedsuch that the center of mass is located at the origin of the coordinatesystem, which simplifies the mathematical notation. With that in mind,the rectangular radial window can be defined as follows:

$\begin{matrix}{{\omega_{r}(d)} = \{ {\begin{matrix}{1,{{d} \leq r}} \\{0,{{d} > r}}\end{matrix};{d = \{ {x,y,z} \}}} } & (7.2)\end{matrix}$

In the one-dimensional case, this degenerates to the familiar recto ( )function, where the width of the rectangle is 2r. In two dimensions, thewindow is circular with radius r. Similarly, in the three-dimensionalcase the window is spherical with radius r. This function, while easy toimplement has several disadvantages. First, it imposes a hard cut-offpoint, giving zero weight to any voxels outside of r. More important,there is no scaling given to any of the interior region. This, ineffect, imposes a hard segmentation based on distance from the center ofmass.

An improvement over the rectangular window would be one that provides agradual attenuation of the voxel intensities as they appear further fromthe COM. Such a window would allow the weighting function to retain themajority of the density information near the center of ROI and much lessnear the periphery. This lead to the choice of a Gaussian weightingfunction.

The one-dimensional Gaussian may be specified in terms of it's standarddeviation, σ. The intensity at a point relative to the mean, μ, isdefined as

$\begin{matrix}{{\omega_{g}(x)} = {\frac{1}{\sigma\sqrt{2\pi}}{{{\mathbb{e}}^{c}( {x - \mu} )}^{2}/2}\;\sigma^{2}}} & (7.3)\end{matrix}$

In two dimensions, the distribution may have a different standarddeviation for each dimension (resulting in an elliptical Gaussian), butas we wish to obtain a symmetric weighting function, we will considerthe circular case whereσ=σ_(x)=σ_(y)The corresponding two-dimensional Gaussian function may be expressed as

$\begin{matrix}{{w_{g}( {x,y} )} = {\frac{1}{2\;\pi\;\sigma^{2}}{\mathbb{e}}^{{{- {\lbrack{({x - \mu})}_{x}^{2}\rbrack}}/2}\;\sigma^{2}}}} & (7.4)\end{matrix}$where the origin is located at (μ_(x),μ_(y)). Similarly, we may define asymmetric (spherical) three-dimensional Gaussian weighting function as

$\begin{matrix}{{w_{g}( {x,y,z} )} = {\frac{1}{{\sigma^{3}( {2\;\pi} )}^{3/2}}{\mathbb{e}}^{{{- {\lbrack{{({x - \mu_{x}})}^{2} + {({y - \mu_{y}})}^{2} + {({z - \mu_{z}})}^{2}}\rbrack}}/2}\sigma^{2}}}} & (7.5)\end{matrix}$Whereas we will set the mean (origin) of the Gaussian weighting functionto be the center of mass of the ROI to be weighted, the choice of anappropriate value for the standard deviation, σ, is less obvious. Theoverall goal of the weighting is to provide a representation of theoriginal ROI with structures at the periphery smoothly attenuated, sothat they contribute less to the overall density distribution. Thus, thechoice of values for σ must be a function of two parameters: the size ofthe ROI, and the degree to which the periphery should be attenuated. Inaddition, the weighting function used for an ROI containing a nodule inone CT scan should be the same as that used to weight the ROI in asubsequent scan. A discussion of the appropriate mechanism for selectionof values for σ will be given in the context of the doubling timeestimation problem in Section 7.3.1.

Moment-Based Analysis—The method of moments may be employed to study thedensity distribution of the weighted ROI. Moments were discussed inSection 6. In this context, we are using densitometric moments asdescriptors of the weighted region containing the nodule. First andforemost is the zeroth-order moment:

$\begin{matrix}{m_{000} = {\sum\limits_{x = 0}^{M - 1}\;{\sum\limits_{y = 0}^{N - 1}\;{\sum\limits_{z = 0}^{L - 1}\;{\upsilon( {x,y,z} )}}}}} & (7.6)\end{matrix}$where v(x, y, z) is a continuous, real-valued function describing the CTattenuation at a given location in 3-space. Recall that when v (x, y, z)is binary-valued, as in a segmented ROI, the zeroth-order momentcorresponds to the volume of the nodule region. In this case, however,it is the sum of the attenuation at each voxel in the ROI. Morespecifically, if we assume that CT attenuation corresponds to objectdensity at a particular location, m₀₀₀ may be thought of as a measure ofregion mass.

In the current application, however, we are analyzing a weighted noduleROI. The goal, then, is to determine a moment-based descriptor that canbe used to characterize the nodule or change in the nodule over time. Wewill, therefore, define the mean density of the weighted nodule ROI as

$\begin{matrix}{{m\; d} = {\frac{m_{000}}{( {M \cdot N \cdot L} )} = \frac{m_{000}}{\mu_{0}}}} & (7.7)\end{matrix}$

This is simply the “mass” of the region divided by the volume of theROI, analogous to physical density measurements. Using this metric, weare able to study the change in size of the nodule over time.

Doubling Time Estimation—Doubling time estimates based on the change insegmented nodule volume were discussed in Section 5. An alternativeapproach to doubling time estimation not requiring explicit segmentationwould be to quantify the change in mean density of the ROI as asurrogate for change in nodule volume.

Consider two scans of an growing object at two times within ROIs ofidentical size. At baseline, the volume of the object is V₁, and in thesecond scan it is V₂. If we consider the case where the object is ofuniform density and that the background is zero-valued (empty), theratio of the mean densities of these regions is exactly equal to theratios of their volumes:

$\begin{matrix}{{{md}_{1} = \frac{V_{1}}{V_{ROI}}},{{md}_{2} = { \frac{V_{2}}{V_{ROI}}\Rightarrow\frac{m\; d_{1}}{m\; d_{2}}  = \frac{V_{1}}{V_{2}}}}} & (7.8)\end{matrix}$

This mathematical relationship suggests that it may be possible to usethe mean density metric as a proxy for nodule volume in computingdoubling time estimates. Of great importance, obviously, is theselection of two corresponding ROIs with identical dimensions. There aretwo additional considerations, however: (a) the nodule may not be ofuniform density; and (b) non-nodule objects may be included in the ROI.

The first consideration supposes that if the nodule changes in densitybut not in “size,” the relationship between mean density and segmentedvolume no longer applies. This is true. However, it may be argued thatan increase in measured nodule density, for equivalent nodule “size,”may be a better indicator of increased numbers of cells in the region.Furthermore, the removal of the intensity-based thresholding step intraditional nodule segmentation allows us to consider the entire nodulevolume in every case, eliminating the complexity and algorithmic (orradiologist) bias in segmentation.

The second consideration, pertaining to the inclusion of non-noduleregions of significant attenuation in the ROI is also valid. However,most of the non-nodule structures in the ROI correspond to vasculatureor small airways. In each case, the size of such structures is likely toremain constant between scans, unless the nodule (as in the case of somemalignant tumors) is responsible for angiogenesis in the region. In thiscase, the additional density measured in the ROI would contribute to amore aggressive estimate of growth, appropriate in such cases.Furthermore, the contribution of structures at the periphery of the ROIwill be limited as the region has been pre-weighted using a Gaussianweighting function. This will help achieve the goal of keeping thedensity distribution of the nodule itself as the chief contributor tothe density distribution of the entire ROI.

Methods—A complete method for nodule doubling time estimation withoutexplicit nodule segmentation has been developed. The goal is to use themean density of the Gaussian-weighted nodule ROI as a proxy for the truenodule volume. In this way, we call compare the mean densities of twoappropriately weighted nodule ROIs taken at two different times anddetermine an estimate of the rate of growth of the nodule. A keyparameter in this procedure is the standard deviation, σ, used in eachof the ROI weighting functions. It is desirable to have the Gaussiantaper toward zero near the periphery of the ROI so that additional,non-nodule structures (e.g. vessels) present in the ROI will have lessof a contribution toward the density calculations. This would suggest arelatively small value of σ might be appropriate. However, we do notwish to over-attenuate the true nodule region. Furthermore, as ourgrowth estimation technique is based on measures of mean density in theregion, we do not wish to choose a value of σ small enough such that themean density metric will become “saturated” upon subsequent nodulegrowth. This might otherwise be possible as we are using the sameweighting function for the ROI at each time. For example, if the firstROI is weighted such that the mean density is quite high (near theeffective maximum of the metric), nodule growth that appears in thesecond ROI may be missed if the radius of the Gaussian is too small.Furthermore, increases in mean attenuation near the center of the nodulein the second ROI may not be accurately measured in this case, as theremight be insufficient dynamic range in the weighted density scale.

For these reasons, an iterative technique is employed to determine thevalue of σ for the Gaussian weighting function such that the meandensity in the first ROI is fixed at a pre-determined point. In thisway, a sufficient amount of dynamic range in the mean density metric ispreserved for measurement of growth in the second ROI Algorithm 7.3describes the method for selection of the appropriate value of σ for agiven standardized mean density, md_(s), and value of ε, the stoppingcriterion.

Algorithm 7.3 (Iterative Sigma Determination) Select an ROI, ROI, and astandardized mean density value, md_(s) σ_(b) = σ_(c) = 0 s = 0.1 md, =0 md_(c) = 0 while done ≠ 1 Generate Gaussian weighting function w_(c),using σ_(c) WROI_(c)(x, y, z) = w_(c)(x, y, z) · (x, y, z)${md}_{c} = {( {\sum\limits_{x = 0}^{M - 1}{\sum\limits_{y = 0}^{N - 1}{\sum\limits_{z = 0}^{L - 1}{{WROI}_{c}( {x,y,z} )}}}} )/( {M \cdot N \cdot L} )}$δ = (md_(s) − md_(c)) if |δ| > 0 σ_(b) = σ_(c) if |δ| < ε done − 1 elseσ_(c) = σ_(c) + s end else if|δ| < ε done − 1 else s = s/2 σ_(c) =σ_(c) + s end end end

At the termination of the algorithm, the value of σ is the appropriatestandard deviation for the Gaussian weighting function, such that itwill produce the desired mean density, md_(s).

Using the mean density of the weighted nodule ROIs as a surrogate forthe nodule volume, it is possible to estimate the growth rate of anodule in sequential CT scans without performing explicit segmentation.Given two 3D nodule ROIs, ROI₁ and ROI₂, acquired at t₁ and t₂,respectively, we may use the first to compute the appropriate weightingfunction for the pair, and subsequently produce the Gaussian-weightedregions of interest, WROI₁ and WROI₂. In each of the Gaussian-weighted3D regions, we then compute the mean voxel densities, md₁ and md₂. Giventhese mean density measures and the time between scans, Δt, we maycompute a region-based doubling time estimate, DT_(R) as

$\begin{matrix}{{D\; T_{R}} = \frac{\ln\;{2 \cdot \Delta}\; t}{{In}( {{md}_{2}/{md}_{1}} )}} & (7.9)\end{matrix}$

This is the analogous expression to the calculation of nodule doublingtime based on segmented volumetric data (Equation 5.24). The entiredoubling time estimation technique is described in Algorithm 7.4.

Algorithm 7.4 (Doubling Time Estimation Without Explicit Segmentation)Select two ROIs, ROI₁, ROI₂, at t₁ and t₂, respectively Align (ROI₁,ROI₂) using Algorithm 7.2 Determine σ based on ROI₁ using Algorithm 7.3Generate Gaussian weighting function w using σ W ROI₁(x, y, z) = w(x, y,z) · ROI₁(x, y, z) W ROI₂(x, y, z) = w(x, y, z) · ROI₂(x, y, z)$ {{md}_{2} = {{( {\sum\limits_{x = 0}^{M - 1}{\sum\limits_{y = 0}^{N - 1}{\sum\limits_{z = 0}^{L - 1}{{WROI}_{1}( {x,y,z} )}}}} )/M} \cdot N \cdot L}} )$$ {{md}_{2} = {{( {\sum\limits_{x = 0}^{M - 1}{\sum\limits_{y = 0}^{N - 1}{\sum\limits_{z = 0}^{L - 1}{{WROI}_{2}( {x,y,z} )}}}} )/M} \cdot N \cdot L}} )$Δt = t₂ − t₁${DT}_{R} = \frac{{ln2} \cdot {\Delta t}}{\ln( {{md}_{2}/{md}_{1}} )}$

It is important to note that this algorithm was designed chiefly tooperate on well-circumscribed and vascularized nodules. Nodules withsufficient proximity or connectivity to the pleural surface (pleuraltail, and juxtapleural nodules) present significant ROT registrationproblems. In particular, both the COM-based and correlation-basedalignment techniques are quite sensitive to ROT selection when asignificant portion of the pleural wall is included.

EXAMPLE 12

The techniques for three-dimensional region analysis without explicitnodule segmentation were compared with those based on segmentation. Thefollowing is a summary of the results.

Doubling Time Estimation—A study was performed using eleven pulmonarynodules for which two high-resolution studies were available. Five weremalignant and six benign, either biopsy-proven or exhibiting no growthover two years. Volumetric doubling times (DT) were computed based onthree-dimensional segmentations of each nodule at each time. Doublingtime estimates (DT_(R)) were also computed using the change in meandensity of the registered ROI, using the techniques described above. Theresults are shown in Table 7.1.

TABLE 7.1 Comparison of in vivo doubling time estimates with and withoutexplicit segmentation Case DT_(R) DT Status 1 193.5 104.4 Malignant 246.70 51.09 Malignant 3 182.8 73.42 Malignant 4 341.4 177.4 Malignant 5224.2 87.04 Malignant 6 4557 825.5 Benign 7 580.2 2026 Benign 8 −2894−1571 Benign 9 1167 395.6 2YNC 10 7162 3335 Benign 11 801.8 845.9 Benign

Comparative statistics were performed to assess the degree ofrelationship between the two types of doubling time estimate. As it wasdifficult to make the assumption both DT and DT_(R) would beGaussian-distributed variables, the standard Pearson correlation wassubstituted with the Spearman rank correlation. A non-parametric measureof correlation, it is based on the difference in rank of correspondingvalues in each variable. If we define R_(i) to be the rank of onevariable in its distribution and S_(i) to be the rank of the secondvariable in its distribution, the Spearman rank-order correlationcoefficient, r_(s) is defined as

$\begin{matrix}{r_{s} = \frac{\sum\limits_{i}{( {R_{i} - \overset{\_}{R}} )( {S_{i} - \overset{\_}{S}} )}}{\sqrt{\sum\limits_{i}( {R_{i} - \overset{\_}{R}} )^{2}}\sqrt{\sum\limits_{i}( {S_{i} - \overset{\_}{S}} )^{2}}}} & (7.10)\end{matrix}$This expression can be simplified if we define the sum squareddifference of ranks, D as

$\begin{matrix}{D = {\sum\limits_{i = i}{N( {R_{i} - S_{i}} )}^{2}}} & (7.11)\end{matrix}$The expression for r_(s) is then simply

$\begin{matrix}{r_{s} = {1 - \frac{6\; D}{N^{3} - N}}} & (7.12)\end{matrix}$The value of r_(s) was computed to show the relationship ofvolumetrically determined doubling times to those determined from changein mean density. In this experiment, r_(s) was equal to 0.9091,P>|r_(s)|=0.0001. This indicates a very high correlation between thevalues of DT and DT_(R). Thus, the estimation of nodule doubling timewithout explicit segmentation may be an effective measure in place ofthe segmentation-based method.

An examination of the DT_(R) values alone as a predictor of nodulemalignancy was also encouraging. The density-estimated doubling timevalues were significantly different for malignant nodules, 197.73±47.1(mean±sem), vs. benign ones, 2860.7±1059.6 (p<0.01, using Wilcoxon ranksums). Again, these techniques for the estimation of nodule growthwithout explicit segmentation appear promising in the differentiation ofbenign from malignant lesions.

Thus, while there have been described what presently believed to be thepreferred embodiment of the invention, those skilled in the art willrealize that changes and modifications may be made thereto withoutdeparting from the spirit of the invention, and it is intended to claimall such changes and modifications as fall within the true scope of theinvention.

1. A method of estimating volumetric doubling time (DT) of an objectwhich changes size, said object found in a mammalian host, comprising:a. obtaining a measurement of a change of a volumetric characteristic ofan object in a mammal over a time period, said object characterized aschanging in size over a time period, said measurement of said change ofsaid volumetric characteristic being obtained from an object at a firsttime (t₁) and at a later second time (t₂); b. relating said change insaid volumetric characteristic found in step (a) to said time period viaan exponential size change model; and c. determining an estimateddoubling time (DT) by comparing volume change over said time periodcomprising: (i) defining a minimum-detectable size dimension β; and (ii)calculating doubling time using volume measurement at t₂ and a volumemeasurement derived from said minimum-detectable size dimension β, andthe time period (Δt_(s)) between t₁ and t₂.
 2. A method as in claim 1wherein said measurement is obtained by two three-dimensional volumetricassessments of an existing object at said first time (t₁) and again atsaid later second time (t₂).
 3. A method as in claim 2 wherein saidexponential size change model is$\frac{V_{2}}{V_{1}} = {\mathbb{e}}^{\lambda\;\Delta\; t}$ wherein V₁and V₂ are first and second volumetric assessments, respectively, and λis an exponential coefficient relating to the doubling time and$\lambda = {\frac{\ln\; 2}{D\; T}.}$
 4. A method as in claim 3 whereinsaid doubling time (DT) is determined as follows:${D\; T} = {\frac{\ln\;{2 \cdot \Delta}\; t}{\ln( {V_{2}/V_{1}} )}.}$5. A method as in claim 1 which further comprises combining aretardation factor (R) with said standard exponential size change model.6. A method as in claim 1 wherein said measurements are obtained fromthree-dimensional images provided from signals bearing three-dimensionaldata produced from scanning said host.
 7. A method as in claim 6 whereinsaid image has been supersampled.
 8. A method as in claim 6 wherein saidimage has been segmented.
 9. A method as in claim 6 wherein said imagehas been supersampled and segmented.
 10. A method as in claim 6 whereinsaid host is human and said object is a growth.
 11. A method as in claim1 and wherein said growth is a pulmonary nodule.
 12. A method as inclaim 11 wherein said pulmonary nodule is not greater than 10 mm indiameter.
 13. A method as in claim 1 wherein said measurement of volumechange is obtained from an object which is not detected at said firsttime (t₁), but which is detected at said second time (t₂).
 14. A methodas in claim 1 which comprises a determination of doubling timeestimation as follows:${D\; T} = {\frac{\ln\;{2 \cdot \Delta}\; t_{s}}{\ln( {V_{2}/V_{\beta}} )}.}$15. A method as in claim 1 wherein said volumetric characteristic ismean density and said change is change in mean density.
 16. A method asin claim 1, wherein said object is or is not detectable at said firsttime (t₁).
 17. An article of manufacture for estimating volumetricdoubling time (DT) of an object which changes size, said object found ina mammalian host, comprising: a machine readable-medium containing oneor more programs which when executed implement the steps of: a.obtaining a measurement of a change of a volumetric characteristic of anobject in a mammal over a time period, said object characterized aschanging in size over a time period, said measurement of said change ofsaid volumetric characteristic being obtained from an object at a firsttime (t₁) and at a later second time (t₂); b. relating said change insaid volumetric characteristic found in step (a) to said time period viaan exponential size change model; and c. determining an estimateddoubling time (DT) by comparing volume change over said time periodcomprising: (i) defining a minimum-detectable size dimension β; and (ii)calculating doubling time using volume measurement at t₂ and a volumemeasurement derived from said minimum-detectable size dimension β, andthe time period (Δt_(s)) between t₁ and t₂.
 18. An article ofmanufacture as in claim 17 wherein said measurement is obtained by twothree-dimensional volumetric assessments of an existing object at afirst time (t₁) and again at later second time (t₂).
 19. An article ofmanufacture as in claim 18 wherein said exponential size change model is$\frac{V_{2}}{V_{1}} = {\mathbb{e}}^{\lambda\;\Delta\; t}$ wherein V₁and V₂ are first and second volumetric assessments, respectively, and λis an exponential coefficient relating to the doubling time and$\lambda = {\frac{\ln\; 2}{D\; T}.}$
 20. An article of manufacture as inclaim 19 wherein said doubling time (DT) is determined as follows:${D\; T} = {\frac{\ln\;{2 \cdot \Delta}\; t}{\ln( {V_{2}/V_{1}} )}.}$21. An article of manufacture as in claim 17 which further comprisescombining a retardation factor (R) with said standard exponential sizechange model.
 22. An article of manufacture as in claim 17 wherein saidmeasurements are obtained from three-dimensional images provided fromsignals bearing three-dimensional data produced from scanning said host.23. An article of manufacture as in claim 22 wherein said image has beensupersampled.
 24. An article of manufacture as in claim 22 wherein saidimage has been segmented.
 25. An article of manufacture as in claim 22wherein said image has been supersampled and segmented.
 26. An articleof manufacture as in claim 17 wherein said host is human and said objectis a growth.
 27. An article of manufacture as in claim 17 and whereinsaid growth is a pulmonary nodule.
 28. An article of manufacture as inclaim 27 wherein said pulmonary nodule is not greater than 10 mm indiameter.
 29. An article of manufacture as in claim 17 wherein saidmeasurement of volume change is obtained from an object which is notdetected at said first time (t₁), but which is detected at said secondtime (t₂).
 30. An article of manufacture as in claim 17 which comprisesa determination of doubling time estimation as follows:${D\; T} = {\frac{\ln\;{2 \cdot \Delta}\; t_{s}}{\ln( {V_{2}/V_{\beta}} )}.}$31. An article of manufacture as in claim 17 wherein said volumetriccharacteristic is mean density and said change is change in meandensity.
 32. An article of manufacture as in claim 17, wherein saidobject is or is not detectable at said first time (t₁).
 33. A system forestimating volumetric doubling time (DT) of an object which changessize, said object found in a mammalian host, comprising: a processorconfigured to: a. obtain a measurement of a change of a volumetriccharacteristic of an object in a mammal over a time period, said objectcharacterized as changing in size over a time period, said measurementof said change of said volumetric characteristic being obtained from anobject at a first time (t₁) and at a later second time (t₂). b. relatesaid change in said volumetric characteristic found in step (a) to saidtime period via an exponential size change model; and c. determine anestimated doubling time (DT) by comparing volume change over said timeperiod comprising: (i) define a minimum-detectable size dimension β; and(ii) calculate doubling time using volume measurement at t₂ and a volumemeasurement derived from said minimum-detectable size dimension β, andthe time period (Δt_(s)) between t₁ and t₂.
 34. A system as in claim 33wherein said processor is further configured to obtain said measurementby two three-dimensional volumetric assessments of an existing object ata first time (t₁) and again at later second time (t₂).
 35. A system asin claim 34 wherein said exponential size change model is$\frac{V_{2}}{V_{1}} = {\mathbb{e}}^{\lambda\;\Delta\; t}$ wherein V₁and V₂ are first and second volumetric assessments, respectively, and λis an exponential coefficient relating to the doubling time and$\lambda = {\frac{\ln\; 2}{DT}.}$
 36. A system as in claim 35 whereinsaid doubling time (DT) is determined as follows:${D\; T} = {\frac{\ln\;{2 \cdot \Delta}\; t}{\ln( {V_{2}/V_{1}} )}.}$37. A system as in claim 33 wherein said processor is further configuredto combine a retardation factor (R) with said standard exponential sizechange model.
 38. A system as in claim 33 wherein said processor isfurther configured to obtain said measurement from three-dimensionalimages provided from signals bearing three-dimensional data producedfrom scanning said host.
 39. A system as in claim 38 wherein said imagehas been supersampled.
 40. A system as in claim 38 wherein said imagehas been segmented.
 41. A system as in claim 38 wherein said image hasbeen supersampled and segmented.
 42. A system as in claim 33 whereinsaid host is human and said object is a growth.
 43. A system as in claim33 wherein said growth is a pulmonary nodule.
 44. A system as in claim43 wherein said pulmonary nodule is not greater than 10 mm in diameter.45. A system as in claim 33 wherein said measurement of volume change isobtained from an object which is not detected at said first time (t₁),but which is detected at said second time (t₂).
 46. A system as in claim33 wherein said processor calculates a determination of doubling timeestimation as follows:${D\; T} = {\frac{\ln\;{2 \cdot \Delta}\; t_{s}}{\ln( {V_{2}/V_{\beta}} )}.}$47. A system as in claim 33 wherein said volumetric characteristic ismean density and said change is change in mean density.
 48. A system asin claim 33, wherein said object is or is not detectable at said sametime (t₁).